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ABSTRACT 

The  effects  of  point  defect  implantation  in  copper  and  tungsten 
crystal  lattice  have  been  studied  by  computer  simulation  techniques 
Vacancies,  interst it ials ,  and  replacement  impurities  have  been 
created  in  the  first  five  layers  of  the  free  (100)  surface  of  these 
crystals.   The  subsequent  binding  energies  of  these  defects  in 
tungsten  were  compared  with  experimental  temperature  dependent  de- 
sorbtion  peaks,  corresponding  to  binding  energies  of  neon  defects 
in  a  tungsten  crystal.   Interstitial  and  replacement  impurity 
positions  in  the  first  three  to  five  layers  were  found  that  seem 
to  correspond  to  the  experimental  data.   Significant  results  were 
also  obtained  which  were  associated  with  general  surface  effects, 
especially  crowdion  migration. 
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I.   INTRODUCTION 

Extensive  research  has  taken  place  in  the  last  decade  in  the 
area  of  computer  simulation  of  radiation  damage  in  crystal  lat- 
tices.  Two  major  areas  of  simulation  have  been  defined.   "Dynamic 
simulation"  suggests  the  firing  of  an  atom  or  ion  against  a  crystal 
and  the  observation  of  the  resulting  many-body  collisions.   Ex- 
amples include  sputtering  simulation  LI, 2, 33,  in  which  atoms  are 
ejected  from  the  surface  of  an  ion-bombarded  crystal;  and  chan- 
neling simulation  l4j ,  in  which  ions  are  fired  down  open  channels 
in  non-close-packed  structures  such  as  body-centered  cubic  and 
diamond  lattices.   "Static  simulation",  on  the  other  hand,  is  con- 
cerned with  equilibrium  positions  and  energies  for  point  defects 
in  crystals  [5,6,7].   This  latter  area  was  the  concern  of  this 
research.   Specifically,  this  simulation  attempted  to  correlate 
equilibrium  potential  energies  of  point  defects  with  experimentally 
determined  binding  energies  of  point  defects  in  Tungsten  [8,9]. 

A.   HISTORICAL  BACKGROUND 
1 .   The  Problem 

Historically,  all  crystal  dynamics  computer  simulation  has 
been  based  on  the  assumption  that  the  complicated  many-body  problem 
can  be  reduced  to  many  two-body  problems.   This  assumption  has  re- 
peatedly been  shown  to  be  a  valid  one  when  employed  incrementally. 
Incremental  calculations  are  necessary  since  a  complete  solution 
in  closed  form  is  impossible.   Small  time  increments  At  approximated 
the  true  time  differential  dt  of  the  impossible  closed  form 
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solution.   Specifically,  the  desired  type,  size,  and  orientati* 
of  crystal  lattice  is  stored  in  the  corn-outer,  appropriate  inter- 
atomic potentials  are  chosen,  and  all  mutual  forces  between  atoms 
are  calculated,  based  on  the  analytic  potential  functions.   Point 
defects  are  introduced,  and  each  atom  is  then  allowed  to  move  in- 
crementally, based  on  these  forces  and  Newtonian  Mechanics  C 5l . 
Through  proper  choice  of  time  increment  duration,  damping  of  forces 
and  velocities  in  each  time  increment,  and  sufficient  repetition 
of  the  procedure,  realistic  results  are  obtained. 
2.   The  Pioneers  and  Their  Contributions 

—     ■■■'■■■■■■  '  '  '■■  ■■      |  '   ■  -^   ■'"      '   -  -  ■'      ■   ■■'-     ■   "     —■■■■■ 

Pioneering   work    in    this    field  began    in    the   early    1960's 
at    the    Brookhaven   National    Laboratory.      Gibson,    Goland,    Milgram 
and   Vineyard    (GGMV)    L 5l    published   the   results    of   extensive   work 
in    both   static   and   dynamic   simulation.      In   static   simulation   they 
determined   equilibrium  positions    for    interstit ials    and   associated 
potential    energies    of   formation.       In    dynamic    simulation    they    in- 
vestigated momentum  propagation    directions    of    energetic    knock-on 
atoms    (focusing),    collision    chains,    and   related   topics.       They   used 
a    central-difference  method   to   obtain    velocities    and  positions    frcm 
calculated   forces.      All   work  was    done   with   copper,    and   results    were 
correlated  with    experimental    data.      Johnson   and   Brown    (JB)    [6 J    did 
extensive   work    in   static    simulation,    again   with    copper.       They   es- 
tablished  that    only   one   stable  position    exists    for    a    single    face- 
centered  cubic    (FCC)    self-interstitial:    the   < 100)    split    inter- 
stitial.     Johnson   [lo]    later   published   further    work    in    this    area, 
with   formation   and   activation    energies    for    various   point    defects. 
Enginsoy.    Vineyard,    and   Englert    (EVE)    [7]    and  Johnson    [llj 


repeated  most    of   the   earlier    calculations    in   GGMV  and  Johnson 
for    the   body-centered   cubic    (BCC)    case,    based   on  QL    iron.      They, 
too,    established  the   existence   of    only    one    stable    interstitial 
position,    a    (lio)    split    interstitial. 

Girifalco   and   Weizer    (GF)    [l2]    calculated  Morse  Potential 

§..   =    D[exp    {-2a(r..-r    )}-    2    exp{-a(r..-r    }] 
ij  r  v    13      o'  r         v    ij      o 

parameters  for  various  metals,  based  on  experimental  values  for 
the  energy  of  vaporization,  the  lattice  constant,  and  compressi- 
bility.  Resulting  elastic  constants  and  equations  of  state  agreed 
satisfactorily  with  experiment.   Girifalco  and  Weizer  L13j  later 
published  results  of  using  these  Morse  parameters  in  simulating 
vacancy  relaxation  dynamics.   Anderman  Ll4l  used  GW's  technique 
of  calculating  Morse  parameters,  but  instead  of  summing  over  an 
entire  crystal,  (GW  calculated  out  to  the  150th  nearest  neighbor) 
Anderman  found  parameters  as  a  result  of  summing  out  to  second, 
third,  and  fourth  nearest  neighbors,  for  use  in  short-range  approxi- 
mations . 

Harrison  [l,2,3]  has  investigated  sputtering  phenomenon  and 
other  surface  effects  with  a  modified  Brookhaven  model,  the  most 
significant  change  being  the  use  of  an  average  force  method  L 1 5 J 
instead  of  the  central  difference  method  in  integrating  the 
equations  of  motion.   (See  Appendix  C. )   He  has  also  calculated 
repulsive  potentials  of  the  Born-Mayer  type  (V.  .  =  exp(A+Br.  .)) 
for  many  combinations  of  atoms  and  ions  based  on  secondary  elec- 
tron emission,  and  Hartree-Fock  atomic  electron  distributions  L 16 J . 


3 .   The  Potential  Function  Problem 

The  most  difficult  problem  encountered  by  computer  sim- 
ulation has  been  the  proper  choice  of  the  potential  function.   No 
simple  analytic  expression,  based  on  either  theory  or  experimental 
data,  has  ever  been  found  that  completely  describes  crystal  dy- 
namics L17J,  although  many  analytic  expressions  are  partially  cor- 
rect.  The  problem  has  been  three- fold:  First,  present  analytic 
expressions  have  narrow  regions  of  validity,  i.e.,  some  correctly 
describe  atomic  behavior  at  equilibrium  distances,  but  fail  at 
shorter  or  greater  distances.   Second,  some  analytic  expressions 
are  limited  because  they  only  apply  to  interactions  between  iden- 
tical atoms.   Third,  the  assumed  functions  have  spherical  symmetry, 
and  are  technically  limited  to  interactions  between  closed  shell 
atoms  or  ions  Ll-J  -   Although  our  assumption  of  a  s^her ic? 13 y  sym- 
metric potential  in  crystals  is  only  approximately  correct,  it  is 
nevertheless  a  very  good  approximation  for  FCC  structures,  and  a 
reasonably  good  approximation  for  BCC  structures.   It  is  grossly 
in  error  when  applied  to  diamond  structures. 

The  atomic  potential,  with  the  familiar  potential  well, 
sharply  repulsive  wall  and  gently  attractive  tail,  varies  greatly 
between  different  pairs  of  atoms;   Since  theory  can  give  only  ap- 
proximate parameters  for  this  complete  potential  function,  experi- 
mental data  have  been  used  extensively  in  the  formulation  of  po- 
tentials.  Other  avenues  have  been  opened  by  computer  simulation. 
Since  potential  well  depths  are  typically  on  the  order  of  a  few 
eV,  the  characteristics  of  the  well  can  be  ignored  in  high  energy 
dynamic  simulation.   Low  energy  dynamic  simulation,  and  even  static 
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simulation,    have  also   been    based  upon    this    approximation  with 
useful    results.       Historically    this    is    how   crystal    simulation 
began.      GGMV  L5l    employed  a   purely   repulsive  potential    of    the 
Born-Mayer     ( BM )    type   and   applied   external    forces    on   all    crystal 
boundaries    to   hold   the    crystal    together.      JB   C6]    used  basically 
the   same   technique.      For    improved   equilibrium   studies   a   potential 
with   a    well    was    necessary    so  GW  [ 12,13]    used  a   Morse  potential    in 
their    simulation.      The  Morse   function,    however,    fails   at    strongly 
repulsive   distances.      To   satisfy    the   need  for    a    more   versatile 
potential,    capable   of   handling    both    high    energy    and   near-equili- 
brium   dynamics,    composite  potentials    were    developed,    which    re- 
semble   BM    or    Bohr 

r  ' 
(V.  .    =  exp(A+Br  .  .)  ) 

functions  at  short  separations,  and  Morse  functions  at  equili- 
brium and  greater  separations.   Specifically,  EVE  L7J  combined 
a  screened  Coulomb  or  Bohr  potential,  a  BM  potential,  and  a  Morse 
potential,  in  the  higher  repulsive,  lower  repulsive,  and  attractive 
regions,  respectively,  of  the  atomic  potential. 

Johnson  [lO,ll]  in  his  later  papers,  used  three  cubic  equations  to 
approximate  the  true  potential.   Anderman  [  l4l  and  Harrison  Ll,2, 
3,4]  have  used  the  BM  repulsive  term  together  with  a  Morse  well 
and  attractive  tail,  smoothly  fit  together  by  a  cubic  equation  in 
the  region  near  their  intersection. 
4.   The  Point  Defect  Problem 

All  early  simulation  was  done  with  homogeneous  systems: 
All  atoms  were  exactly  the  same,  limiting  high  energy  dynamic 
simulation  to  bombardment  by  atoms  identical  to  the  lattice  atoms, 
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and   static    simulation    to   consideration    of   only    the    vacancy   and 
self-interstitial    cases.      This    limitation  was    forced  by    the  po- 
tential   functions,    because  parameters    for    the  Morse   function   were 
based   on    experimental    data    for    homogeneous   media    [12],      The 
methods    used  could   not   yield  parameters    for    different-atom  pairs. 
BM  parameters,    however,    are    obtainable    for    different-atom  pairs, 
by   methods    such   as    the   Hartree-Fock   method   Cl6J    mentioned  pre- 
viously . 

In    spite  of   these    limitations,    dynamic    computer    simulation 
of  bombardment    by   foreign   atoms,    or    static    simulation    of   foreign 
interst it ials ,    can   be    done   by   two   alternate   methods.       First, 
Harrison   [4l    has    neglected   the  attractive    interactions    and   has 
done   foreign  particle    dynamics    using   repulsion    only.      Alternately, 
Jchnscn    111]    has    derived   a    cubic    equation    for    a    complete   DOtential 
with   a  potential    well,    based   on    limited   experimental    data    on    car- 
bon   defects    in    iron.       However,    experimental    substantiation    for    a 
foreign-particle  potential    well    is    much   more    difficult    than    for 
an    identical    atom  potential   well. 

B.       THE    EXPERIMENT 

The    experimental    data    which   this    simulation   proposed  to    explain, 

were  published  by   Kornelsen   and   Sinka    (KS)    L8].      They   have  bom- 

+  +  +  + 

barded  a    clean    (100)    tungsten    surface   with   Ne    ,    Ar     ,    Kr     ,    and   Xe    , 

in    the    energy    range   of   40   eV    to   5   keV.      The   subsequent    "damaged" 

crystal    was    heated  at   a    constant    rate,    and   gas    desorbtion    rates 

were  measured.       Instead   of   a    constant    desorbtion    of    ions,    various 

distinct   peaks    were   found,    categorized    into   two  basic   types:    a 
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single  large  peak  at  1800  k,  the  same  for  all  four  ions;  and  four 
or  five  smaller  peaks  in  the  400  K  to  1650  K  range,  which  were 
not  in  the  same  position  for  all  four  ions.   These  latter  peaks 
were  postulated  to  correspond  to  binding  energies  of  various  point 
defects  in  the  first  few  layers  of  the  tungsten  crystal.   (See 
Figure  4. ) 

This  simulation  used  Harrison's  assumption  of  a  repulsive 
potential  only,  for  interactions  between  a  foreign  point  defect 
and  other  atoms  in  the  lattice.   When  investigating  neon  defects 
in  tungsten,  all  tungsten  lattice  interactions  were  based  on  com- 
posite Morse  and  repulsive  Born-Meyer  potentials,  and  all  neon- 
tungsten  (Ne-W)  interactions  were  based  on  a  purely  repulsive  BM 
potential . 
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II.   OBJECTIVE 


The  long-range  objective  of  this  simulation  was  to  correlate 
simulated  and  experimental  binding  energies  of  neon  point  defects 
in  tungsten.   Since  the  assumption  that  all  Ne-W  interactions  are 
purely  repulsive  was  not  realistic,  the  degree  to  which  subsequent 
simulation  results  are  valid  must  be  based  on  a  known  standard.   If 
the  simulated  binding  energies  are  not  correct,  a  valid  correction 
factor  can  be  applied,  if  derivable  from  the  known  standard.   One 
standard  which  proved  to  yield  this  information  was  the  tungsten- 
tungsten  (W-W)  interaction.   A  tungsten  point  defect  could  be 
treated  as  an  atom  of  the  lattice,  and  given  an  interatomic  po- 
tential identical  to  all  other  lattice  atoms,  the  composite  Morse 
and  BM  potential.   This  was  Method  1.   A  tungsten  defect  could 
also  be  treated  as  foreign,  and  allowed  to  interact  with  other 
lattice  atoms  with  a  repulsive  potential  only  (Method  2).   If  a 
specific  tungsten  point  defect  is  treated  by  both  of  these  methods, 
an  empirical  relationship  between  the  repulsive  potential  assumpt- 
ion and  the  "true"  potential  for  W-W  interactions  is  obtained. 

The  objectives  of  this  research  were  fourfold: 

1.  Demonstrate  that  the  two  methods  of  treating  tungsten  point 
defect  in  a  tungsten  lattice  yield  basically  the  same  physical 
results,  and  agree  with  published  results  [7,ll]  concerning 
split  interstitial  positions. 

2.  Develop  a  general  empirical  relationship  between  the  binding 
energies  derived  by  the  two  methods. 
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3.  Obtain  values  for  binding  energies  of  neon  defects  in  all 
possible  positions  in  a  tungsten  surface.   Transform  these 
values  to  more  realistic  ones  using  the  empirical  relationship 
derived  in  2. 

4.  Compare  these  results  with  KS's  experimental  data. 
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III.       THE    MODEL 


A.       THE    CRYSTAL 

The   model   used    in    this    research    is    the  Gay-Harrison   [lQ]    model, 
with   modifications    by   Levy    [20],    Johnson    [2l],    Effron    [22],    and 
Moore  L23].      Abbreviations    in    brackets    refer    to   computer    program 
names    for    the   variable    in    question. 

Both    copper    and  tungsten    crystals    were   simulated.      Copper    was 
simulated   only    to  provide   an    interface   between    this    research  and 
published  simulation    results.       Copper    forms    a    face-centered  cubic 
crystal    with   an    experimentally    determined   lattice   constant,    (LC) 
or    cube   edge    distance   of   3«6l5A.      The    lattice   unit    (LU),    defined 
as   %LC.    is    1.8075A:    and   the   nearest    neighbor    distance,    as    in   all 
FCC  structures,    is     \T~2LU .      Tungsten    forms    a    body-centered   cubic 
crystal,    with   a   LC  of   3-16A,    a    LU   of    1 . 58A ,    and   a   nearest    neighbor 
distance,    peculiar    to   all    BCC   structures,    of  \|3lU.      All    distances 
in    the  program   are  measured   in   LU.      The  program   could   construct 
(100),    (110),    and    (111)    orientations    of   face-centered  and   body- 
centered   cubic    structures.      The   copper    crystal    size   was    8x8x8, 
and   contained   256   atoms    for    the    (100)    orientation. 

The  major    portion    of    the   simulation   was    done   on    the    (100) 
orientation    of   tungsten,    corresponding    to  KS's    experimental    work. 
This    tungsten    crystal    size    for  Neon   point    defects    was    10   x    10   x    10, 
and   contained  250   atoms.      Some  W-W  simulation   was    done    on   a 
l4   x    14   x   l4    crystal;    the   reasons    are   explained    in    RESULTS .      The 
bottom   two   layers    of   the    lattice  were   not   allowed   to   move,    al- 
though   they   had  potential    energy,    and    exerted   force   on   all   atoms 
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in    the    crystal.      The   other    eight    layers    were   completely    free   to 
move,    and  were    included    in    the   dynamic    calculations    of   each    time- 
step. 

The   surface    layer    (Y   =   0)    and   the   second    layer    (Y    =    1)    were 
moved  forward,    simulating   actual    surface   relaxation    in    the   crystal. 
This   relaxation   was    calculated   by  Moore   [24]    using    simulation 
techniques,    and   tested  against   previous    results    by    Burton   and 
Jura   L"25]-       Definitions    and  use    of  mobile    layers,    relaxation,    etc., 
were  analagous    in    the    copper    model,    as    were  all    other    aspects    of 
the  model    to  be   described   in    this    chapter. 

B.      THE  POTENTIALS 

1.    The  W-W   Composite  Potential 

The  attractive  potential   used   was    the  Morse  potential,    with 
tungsten  parameters    calculated  by   GW   L 12J .      The    interaction    energy 
0.  .,    of   a   pair    of  particles    i   and    j    is: 

0.  .    =    D[exp{-2a(r  .  .-r    )}    -    2    exp[-a(  r  .  .-r     )}]  (1) 

ij  r  v     ij      o'  r         v     ij       o'  v 

where    D  [dC0n]    is    the    dissociation    energy    of   the  pair,    r    [re]    is 
the    equilibrium   separation,    r.  .[DISTj    is    the  actual    separation, 
and  Ci   [ALPHA]    is    a    constant. 

The   repulsive  potential    is    of    the    BM   type,    with   Harrison's 
Hartree-Fock  parameters.      The    interaction    energy   V. .,    is 

V.  .    =    exp    (A+Br  .  .)  (2) 

where  B  [exb]  is  always  negative,  A  [exa]  is  always  positive,  and 
r.  .  [diSTJ  is  the  actual  separation.   The  constants  [exa]  and 
[exb]  are  peculiar  to  the  W-W  interaction. 
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The  ranges  of  the  W-W  composite  potentials  were  as  follows: 
the  BM  repulsive  potential  operated  from  0  to  1 .  5A;  and  the  Morse 
potential  from  2$  [ROES]  to  5.38A*  Croec] .   In  LU ,  the  dimension 
in  which  all  calculations  were  done,  these  constants  were  -9494, 
1.2658,  and  3.4000  LU.   CROEC]  was  chosen  to  include  interactions 
out  to  the  fourth  nearest  neighbor  (NN4)  at  \[Tl      LU  =  3-317  LU  but 
not  NN5  interactions,  at  \|T2  LU  =  3-464  LU.   Note,  however,  that 
slight  displacements  of  NN5's  might  alloiv  their  inclusion  in  po- 
tential and  force  calculations.   The  gap  between  1 . 5A  and  2A  was 
filled  with  a  cubic  function,  which  matched  to  the  other  po- 
tentials and  slopes  at  [ROEA]  and  [r0Eb]. 
2 .   Purely  Repulsive  Potentials 

For  foreign  point  defect  interactions,  i.e.,  Ne-W,  or 
W-W  Method  2.  a  repulsive  potential  only  was  used.   The  potential 
was  again  a  BM ,  with  the  constants  labeled  [PEXAJ  and  Cpexb] .   For 
the  W-W,  Method  2  interaction,  LpEXA]  =  [exa'J  and  [>EXB]  =  [exb] . 
Ranges  for  foreign  point  defect  interactions,  however,  were  dif- 
ferent, and  the  potential  itself  was  modified  at  the  cutoff  point. 
Whereas  the  BM  part  of  the  composite  potential  extended  out  to 
about  .95  LU,  the  modified  BM  potential  used  for  foreign  defects 
was  allowed  to  extend  to  \|3  LU,  corresponding  to  the  NNl  distance 
[ROE].   Cutting  the  potential  off  at  [roe]  left  a  step  of  about 
.05  eV  for  Ne-W  (.2  eV  for  W-W)  at  the  NNl  equilibrium  position. 
Since  neither  discontinuities  nor  repulsive  potentials  were  de- 
sired at  this  equilibrium  position,  we  "eroded"  L15J  the  potential 
by  subtracting  V([R0E]),  or  about  .05  eV  from  V(r. .)  for 
r.  .  <  [roe].   Calculated  forces,  based  on  these  eroded  potentials 
must  be  modified  also,  but  for  a  different  reason.   It  is  possible 


to  conceive  of  a  case  where  an  atom  is  further  away  from  the 
defect  than  [roe]  at  the  beginning  of  a  timestep,  but  closer  than 
[ROE]  at  the  end.   The  force  has  essentially  "turned  on"  in  the 
middle  of  the  timestep.   The  modification  gives  such  an  atom  a 
force  which  is  less  than  the  final  force  by  approximately  a  fac- 
tor proportional  to  the  ratio  of  distance  traveled  outside  [roe! 
to  the  total  distance  traveled  during  the  timestep  Cl5l.   (See 
Appendix  B.  ) 

C.   THE  TIMESTEP 

Motion  caused  by  these  forces  must  be  found  by  an  approximate 
numerical  method  of  time  integration.   As  in  previous  work  with 
this  model,  the  average  force  method  L 1 5J  was  used.   In  this 
method,  all  mutual  forces  were  calculated  in  subroutine  STEP. 
Based  on  these  forces,  new  temporary  velocities  and  positions 
were  found.   Forces  were  again  calculated,  based  on  the  temporary 
positions.   The  final  positions  were  then  calculated  from  the 
average  of  these  two  force  determinations.   All  velocities  were 
then  either  zeroed  or  halved  as  a  damping  method.   This  consti- 
tuted one  timestep. 

For  this  average  force  method  to  work  properly,  the  Atforl 
which  approximates  dt  in  the  integration  must  be  kept  small.   Too 
small  a  value  for  Ldt],  however,  would  result  in  excessive  com- 
puter time.   The  choice  of  LDTj  was  also  complicated  by  the  fact 
that  Cdt]  must  be  kept  much  smaller  earlier  in  the  program,  when 
velocities,  forces,  and  energies  are  large,  but  can  be  allowed  to 
grow  larger  as  the  simulation  approaches  equilibrium.   For  this 
reason,  at  the  end  of  each  timestep,  a  new  LDT]  was  calculated 
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for  use  in  the  next  timestep.   The  parameter  chosen  to  control 
Cdt]  was  [DTl],the  distance,  measured  in  LU,  which  the  most 
energetic  atom  was  allowed  to  move  before  starting  a  new  timestep. 
LCTlJ  has  varied  between  .001  and  .02,  depending  on  such  con- 
ditions as  original  position  of  the  point  defect,  relative  masses 
of  atoms,  etc.   In  general,  L  DTll  must  be  kept  very  small  when 
high  velocities  are  expected,  and  can  be  increased  when  all  motion 
is  expected  to  be  "sluggish".   In  actual  practice  [dt]  and  [DTl] 
are  related  to  both  the  velocity  of  each  particle  and  the  force 
on  each  particle.   To  insure  that  no  particle  traveled  more  than 
CdTI]  we  ensured  that  LDX]  was  small  enough  so  that  neither  the 
velocity  of  the  most  energetic  atom  nor  the  force  on  the  most 
.stressed  atom  would  result  in  motion  greater  than  LDTlJ.   (See 
Appendix  C. ) 

D.   FOREIGN  INTER ST IT I ALS 

1 .   Unequal  Mass  Implications 

Many  changes  in  the  program  were  necessary  when  a  foreign 
defect  was  included  in  the  lattice.   These  changes  were  especially 
necessary  when  the  defect  was  much  lighter  than  the  lattice  atoms; 
i.e.,  neon  in  a  tungsten  lattice.   First,  in  the  average  force 
calculations,  a  separate  section  had  to  be  added  for  the  calcu- 
lations for  the  primary,  or  "bullet",  based  on  the  bullet  mass 
CbMAS] .   Second,  the  potential  energy  between  two  unequal  mass 
atoms  was  split  in  proportion  to  their  reduced  masses  (see 
Appendix  D)  .   Third,  the  section  that  determined  the  new  timestep 
duration  was  originally  based  en  the  lattice  atom  mass.   Since  the 
light  interstitial  is  usually  the  most  energetic  or  most  stressed 
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atom,  very  erratic  behavior  was  observed  until  the  timestep  du- 
ration calculations  were  revised  to  handle  two  different  masses 
(see  Appendix  C) .   Finally,  a  significant  mass  difference  between 
defect  and  lattice  atom  required  a  reduction  in  TdTI] .   For  the 
neon-tungsten  simulation,  a  [DTl]  of  .5%   was  used. 
2.   Ionization  State  and  Repulsive  Potentials 

The  major  portion  of  the  Ne-W  work  was  done  with  the 
assumption  that  tungsten  was  in  a  +6  state,  and  neon  was  neutral 
in  the  lattice.   Experimentally  KS  fired  neon  in  a  +1  state  into 
tungsten,  but  once  emplanted,  the  ionization  of  neon  was  unknown. 
All  combinations  of  W  ,  W   and  W    with  Ne   and  N  '  were  sub- 
jected to  Hartree-Fock  analysis.   (See  Figure  5-)   Only  Ne  -W 
interacted  in  an  approximately  exponential  manner  and  could  there- 
fore possess  realistic  BM  parameters.   Attempts  to  linearize 

o   +1    ,    o   o 
Ne  -W   and  Ne  -W   were  made,  and  subsequent  BM  parameters  were 

determined.   The  results  of  such  changes  did  not  significantly 

influence  the  results  of  this  investigation. 

E.   RUNNING  TIME 

The  following  factors  effected  the  problem  running  time:  range 
of  potential,  size  of  crystal,  depth  of  mobile  layers,  and  degree 
of  damping.   First,  the  range  of  the  potential  was  picked  to 
include  at  least  NN4  interactions.   The  range  used  for  the  tungsten 
simulation  was  3-4  LU,  which  includes  interactions  out  to  NN4 •   The 
error  made  by  neglecting  NN5  interactions  was  only  3%  in  the 
binding  energy  of  an  interstitial,  but  the  omission  of  NN5  inter- 
actions cut  running  time  almost  10%.   Second,  the  size  of  crystal 
and  depth  of  mobile  layers  were  picked  as  small  as  possible,  for 
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reduced  running  time,  but  were  at  least  large  enough  to  completely 
contain  the  potential  range.   Third,  the  half  velocity  method  of 
damping  was  used  whenever  possible.   In  general,  when  velocities 
were  zeroed  at  the  end  of  each  timestep,  a  timestep  took  about  ten 
seconds,  and  equilibrium  was  reached  in  about  300  timesteps.   When 
velocities  were  halved,  each  timestep  again  took  about  ten  seconds, 
but  equilibrium  was  reached  in  about  150  timesteps. 

It  was  originally  expected  that  increasing  LdTIJ  would  decrease 
running  time.   Most  W-W  simulation  was  done  with  LDTlJ  =  2%;  i.e., 
the  most  energetic  or  stressed  atom  could  travel  .02  LU  before  the 
damping  of  velocities  and  the  starting  of  a  new  timestep.   When 
fDTl]  was  increased,  the  atoms  moved  more  erratically  toward  equi- 
librium, and  vibrated  there,  but  did  not  achieve  equilibrium 
significantly  scorer .   (See  Figure  6.) 

F.   SUMMARY 

In  summary,  the  steps  of  the  program  are  outlined: 

1.  Variables  are  initialized,  constants  established,  and  input 
data  read  in. 

2.  Scaling  factors  and  time  saving  multipliers  are  calculated. 

3.  Morse  and  BM  potential  functions  are  calculated  based  on 
input  data.   Subsequent  forces,  based  on  derivatives  of  these 
functions  are  calculated.   Potential  erosion  and  force  modifi- 
cations are  performed. 

4.  Potential  cutoff's  [R0Ea] ,  CROEbI ,  [rOECJ  are  established 
and  the  smooth  fitting  cubic  equation  is  placed  in  the  gap. 

5.  The  desired  crystal  type,  size,  and  orientation  is  built, 
and  the  point  defect  positioned  (see  Appendix  A). 
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6.  Mutual  potential  energies  of  all  atoms  in  the  crystal  are 
calculated.   Local  potential  energy  is  calculated.   (See 
Appendix  D. ) 

7-  All  initial  positions  and  potential  energies  are  printed, 
along  with  total  potential  and  total  kinetic  energy,  local 
potential  energy,  and  the  change  in  local  potential  energy. 

8.  The  first  timestep  is  started,  with  an  arbitrary  running 

-14 

time  of  10     sec.   Velocities  and  positions  are  calculated 
by  the  average  force  method,  and  the  maximum  velocity  L EM AX J 
and  maximum  force  IFMAXJ  are  found. 

9.  A  new  [or],  based  on  [EMAX]  ,  CfMAX']  ,  and  [DTI]  is  calculated 
for  use  in  the  next  timestep.   (See  Appendix  C.) 

10.  All  velocities  are  zeroed  or  halved  as  an  energy  damping 
method;  and  the  process  (8.  to  10.)  is  repeated. 

11.  At  selected  timesteps,  all  changes  in  position  ( L DX J ,  CdyJ , 
CdzJ),  velocities  ([vx],  L w] ,  LvxJ),  and  kinetic,  potential, 
and  total  energies  L'pKe]  ,  [pPe]  ,  [pTe]  for  each  atom  in  the 
crystal  are  printed. 

12.  The  program  is  ended  after  a  pre-selected  timestep,  with  a 
final  printout  of  position  and  potential  energy  of  each  atom, 
as  in  Step  7. 
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IV.   RESULTS 

KS's  experimental  data  indicated  that  four  or  five  interstitial 
positions  in  the  first  few  layers  of  a  tungsten  lattice  could  be 
found  that  would  result  in  different  binding  energies.   It  was  soon 
found  that  many  parameters  in  the  program  could  effect  the  results, 
and  so  a  systematic  attempt  to  isolate  the  effects  of  each  indivi- 
dual parameter  was  undertaken. 

A.   THE  CRYSTAL 

As  explained  in  Appendix  A,  the  Y  =  0  plane  was  the  crystal 
surface;  the  Y  =  1  plane  was  the  first  layer  beneath  the  surface, 
etc.   Each  atom  in  each  layer  was  then  designated  by  appropriate 
x  and  z  coordinates.   This  consxruction  was  independent  of  type 
of  lattice;  i.e.,  the  surface  layer  in  either  BCC  (100)  or  FCC 
(100)  was  Y  =  0,  etc.   An  interstitial  that  escaped  the  lattice 
normal  to  the  surface  travelled  in  a  (Oio)  direction,  and  an  ion 
that  escaped  normal  to  a  side  travelled  in  either  a  (l00>  or  (001> 
direction . 

The  dimensions  of  the  lattice  had  a  great  bearing  on  the  re- 
sults.  In  general,  the  larger  the  crystal,  the  more  realistic  the 
results,  but  increased  computer  time  prevented  the  use  of  a  size 
bigger  than  absolutely  necessary.   All  point  defects  were  placed 
as  close  to  the  center  of  each  plane  as  possible,  and  the  x  and  z 
dimensions  of  the  lattice  Cix]  and  [ iz]  were  chosen  to  completely 
enclose  a  circle  of  radius  [rOEC]  from  the  point  defect.   Tn  this 
way  any  point  defect  in  the  center  of  the  lattice  would  not  feel 
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the  effect  of  the  sides  of  the  lattice,  especially  unequal  numbers 
of  atoms  in  all  directions.   In  the  y  direction,  the  lattice  was 
again  built  deep  enough  to  completely  contain  the  radius  of  the 
potential  of  a  point  defect  placed  at  the  center  of  the  lattice. 
To  simulate  the  effect  of  an  infinitely  deep  lattice,  the  bottom 
two  layers  of  the  crystal  were  held  immobile,  but  still  allowed 
to  interact  with  all  mobile  atoms  above  them.   The  tungsten  crystal 
size  used  most  often  was  a  10  x  10  x  10  cube  with  the  bottom 
10  x  2  x  10  volume  held  rigid. 

Often  erratic  behavior  in  the  simulation  could  be  eliminated 
by  simply  increasing  the  crystal  size.  This  was  especially  true 
for  the  problem  of  crowdion  migration. 

B.   CROWDION  MIGRATION 

Crowdion  migration  is  a  chain  reaction  of  single  lattice  site 
jumps  initiated  by  interstitial  implantation.   If  the  chain  re- 
action ends  by  pushing  the  surplus  atom  into  an  already  existing 
vacancy,  the  interstitial-vacancy  pair  is  called  a  Frenkel  pair. 
A  Frenkel  pair  can  also  be  created  dynamically  by  moving  an  atom 
from  its  lattice  site  to  a  nearby  interstitial  position,  from 
where  it  can  cause  migration  back  to  the  vacancy.   If  the  mi- 
gration cannot  find  a  vacancy,  and  travels  all  the  way  to  the 
surface,  the  surplus  atom  forms  a  "stub".   Normally,  migration 
is  always  in  a  closed  packed  direction;  i.e.,  in  the  (ill)  di- 
rection in  BCC.   It  was  discovered,  however,  that  this  rule  was 
modified  near  a  surface,  since  an  imbalance  of  forces  in  the  di- 
rection normal  to  the  surface  automatically  pushed  a  crowdion  in 
that  normal  direction  into  a  stub  position.   In  the  tungsten 
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lattice,    for    instance,    a    tungsten    interstitial    that    did   not 
initiate   crowdion   migration   would   reach   equilibrium    in   a    (llO> 
split    interstitial   position,    as   previously    found  by   EVE   C7] •      A 
tungsten    interstitial    that    did    initiate   crowdion   migration   would 
sometimes    migrate    in   a    (ill)    direction   because    of   closed-packedness , 
or    sometimes    in   a    (lOO)    direction    if    implanted  near    a    (100)    sur- 
face.     Crowdion   migration  was    never    found    in   a    (llO)    direction, 
since   this    is    the   least    closed-packed  of  these    three   directions. 
(Hence   the    tendency    toward   split-inter st it ials    in   this 
direction ) . 

Crowdion  migration   was    a    very   common   process   near    a    lattice 
surface.       It    was    found,    however,    that    varying    the   choice   of  atoms, 
the  range   of    the  potential,    and    the   rate    of   energy    damping   could 

already   mentioned   was    the   fact    that    increased  crystal    size   re- 
reduced    crowdion   migration.      Particular    attention   was   paid   to   the 
proper    choice   of    values    for    these  parameters,    in    order    to   cor- 
rectly   determine    whether    or    not    crowdion   migration    actually    existed. 
This    question   of   crowdion   migration   was    especially    critical    in    this 
simulation,    since   the  binding    energy    of  a   particular    atom    is    a    di- 
rect   function   of    the    nearness    of    its    neighbors.      An   atom    in   a 
split    interstitial   position    feels    a   more   repulsive  potential    than 
an    interstitial    that    has    initiated  crowdion   migration   and 
'stolen'    a    lattice    site  and    has    thus    reformed  the   original   perfect 
lattice   with   every    atom    in   a    normal    lattice   site. 
1 .       Choice    of   Atom 

Some    elements    tended   to    initiate   crowdion   migration   more 
than    others.      This   applied   to   both    choice   of    lattice   atom,    and 
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choice  of  interstitial  atom.   For  instance,  crowdion  migration 
was  much  more  common  in  tungsten  than  in  copper.   This  was  due 
to  both  the  size  of  the  tungsten  atom  and  the  nature  of  the  crystal 
Also  it  was  found  that  different  element  point  defects  in  the  same 
lattice  produced  varying  degrees  of  crowdion  migration.   A  neon 
interstitial  is  so  small  and  light  that  it  never  initiated  crowd- 
ion migration,  even  when  only  one  layer  separated  it  from  the  sur- 
face.  An  argon  interstitial  initiated  crowdion  migration  at  the 
surface,  but  not  deeper  in  the  lattice.   A  tungsten  interstitial 
always  initiated  crowdion  migration,  unless  placed  at  the  center 
of  a  huge  l4  x  l4  x  14  tungsten  lattice.   This  is  an  example  of 
increasing  crystal  size  to  prevent  crowdion  migration.   In 
general  it  can  be  stated:  the  more  massive  the  interstitial,  the 
more  probable  crowdion  migration. 
2 .   Range  of  Potential 

In  surface  simulation,  the  range  of  the  potential  was  a 
critical  factor,  since  it  determined  whether  or  not  an  atom  could 
"see"  the  surface.   Because  copper  has  been  the  standard  element 
for  lattice  simulations,  many  versions  of  a  copper  potential  with 
various  ranges,  have  been  determined.   As  mentioned  previously, 
GW  [l2]  calculated  a  Morse  potential  for  copper  that  effectively 
had  an  infinite  range  (150th  nearest  neighbor).   If  this  potential 
is  truncated  at  very  close  ranges,  i.e.,  NNl  or  NN2 ,  the  potential 
is  seriously  underestimated.   This  under  estimate  rapidly  dimin- 
ishes as  the  truncation  range  increases.   Since  GW  parameters  for 
the  Morse  potential  could  not  be  used  for  NN2  interactions, 
Anderman  Cl4]  calculated  parameters  for  a  Morse  copper  potential 
that  would  approximate  GW ' s  results  in  simulations  truncated  after 
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NN2 .   He  did  this  by  deepening  and  broadening  the  well.   Although 
Anderman  parameters  and  GW  parameters  led  to  very  similar  results 
in  an  infinite  lattice,  they  led  to  quite  different  results  in  this 
simulation.   In  general,  if  the  range  of  a  point  defect  potential 
function  overlapped  a  surface,  crowdion  migration  would  take  place 
toward  that  surface,  because  of  an  imbalance  of  forces  in  the  nor- 
mal direction.   The  effect  was  a  little  more  complex  than  this  be- 
cause of  surface  relaxation:  if  the  range  of  the  point  defect  po- 
tential function  overlapped  a  relaxed  surface  layer,  the  slight 
force  imbalance  would  again  result  in  crowdion  migration.   Accord- 
ing to  GGMV,  "the  machine  calculation  showed  that  this  atom  rapidly 
moved... in  a  direction  determined  by  minor  asymmetries  in  the 
starting  conditions..."  .   In  this  copper  simulation,  an  inter- 
stitial placed  in  the  forth  layer  with  a  GW  potential  range  of 
3-1  LU  (NN4)  caused  complete  crowdion  migration,  resulting  in  a 
copper  stub  on  the  surface.   An  identical  run  with  Anderman  para- 
meters for  an  NN2  potential  to  a  range  of  2.4  LU  resulted  in  a 
( 100)  split  interstitial  with  minor,  damped  migration  to  the  sur- 
face.  Instead  of  a  stub  copper  atom  as  before,  four  copper  atoms 
in  the  surface  layer  bulged  about  .4  LU. 

Another  example  of  a  short  range  potential  which  demon- 
strated this  lack  of  ability  to  initiate  crowdion  migration  was 
the  repulsive  foreign  defect  potential,  with  a  range  of  \[~3~  LU. 
In  copper,  this  potential  quickly  led  to  split  interstitial 
positions  and  no  crowdion  migration  for  all  copper  interst it ials 
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except  those  placed  in  the  first  two  layers.   In  tungsten,  even 
this  short  range  potential  could  not  retard  crowdion  migration  in 
the  10  x  10  x  10  lattice.   Only  in  the  center  of  a  14  x  14  x  l4  was 
a  tungsten  split  interstitial  stable.   This  stability  applied  only 
to  the  short  range  repulsive  potential:  a  repeat  run  using  the 
standard  composite  potential  with  a  range  of  3.4  LU  initiated 
crowdion  migration.   This  increased  range  enabled  the  interstitial 
to  find  minor  asymmetries  in  even  a  l4  x  l4  x  14  lattice. 
3 •   Energy  Damping 

Energy  damping  was  accomplished  in  this  simulation  by 
reducing  each  atom's  velocity  at  the  end  of  each  timestep.   Two 
methods  were  used:  at  first,  each  velocity  component  of  every 
atom  in  the  crystal  was  zeroed  at  the  end  of  each  timestep.   Later, 
the  halving  of  each  velocity  component  at  the  end  of  each  time- 
step  was  employed  to  save  computer  time.   In  a  tungsten  lattice 
the  results  of  both  methods  were  the  same;  all  final  positions  and 
binding  energies  were  identical.   Neither  method  prevented  crowd- 
ion migration.   In  copper,  these  two  methods  led  to  slight ly ' dif- 
ferent results.   Although  the  final  position  and  binding  energy 
of  an  interstitial  was  almost  identical,  and  although  crowdion 
migration  was  initiated  in  both  cases  (GW's  parameters  and  a 
3.1  LU  range  were  used),  the  zeroed  velocity  method  had  damped  the 
migration  significantly  by  the  time  it  reached  the  surface,  where- 
as the  halved  velocity  method  caused  a  complete,  undamped  mi- 
gration to  the  surface. 
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C.   INTERSTITIAL  IMPLANTATION 

All  inter stitials  were  placed  in  the  obvious  holes  in  a  hard 
sphere,  close-packed  lattice  model.   (See  Figure  7.)   Every 
"hole"  in  the  tungsten  lattice  had  exactly  the  same  geometry; 
i.e.,  two  neighbors  1  LU  away;  four  neighbors  2  LU  away,  four 
neighbors  3  LU  away,  etc.   The  only  factors  which  differentiated 
between  these  identical  holes  and  thus  led  to  different  binding 
energies  were:  layer  number,  or  lattice  depth,  and  open  channel 
direction.   An  interstitial  in  the  third  layer  was  more  tightly 
bound  than  one  in  layer  two,  etc.   Also,  an  atom  in  a  given  layer 
could  be  placed  in  two  types  of  holes:  one  in  which  the  inter- 
stitial was  in  the  BCC  (010)  open  channel  direction,  in  which  case 
the  interstitial  could  "see"  the  surface;  and  one  in  which  the 
interstitial  was  in  the  BCC  ( 100)  or  (OOl)  open  channel,  in  which 
the  interstitial  could  not  "see"  the  surface.   Note,  however,  that 
if  an  atom  could  not  "see"  the  surface,  then  there  was  no  difference 
between  these  two  positions,  since  both  have  two  neighbors  1  LU 
away,  four  neighbors  2  LU  away,  etc.   Note  also  that  even  if  these 
two  sites  are  identical  and  possess  exactly  the  same  binding 
energies,  the  difference  might  still  show  up  in  diffusion  pro- 
bababilities :  an  interstitial  in  a  (010)  open  channel  in  the 
second  layer  must  only  move  two  lattice  units  to  escape  the  cry- 
stal.  An  interstitial  in  a  (lOO)  or  (001)  open  channel  must  move 

1  LU  in  either  the  x  or  z  direction  into  an  open  channel,  and  then 

2  LU  to  escape;  i.e.,  it  must  move  like  a  knight  in  chess.   This 
extra  step  might  lead  to  a  different  diffusion  probability. 
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D.   THE  TUNGSTEN  LATTICE  SELF  DEFECT 

The  tungsten  self  inter stitials  and  self  replacement  defects 
were  the  chosen  standard  for  this  analysis,  as  explained  in  the 
in  the  OBJECTIVE.   The  tungsten  defects  could  be  treated  as 
lattice  atoms  and  allowed  to  interact  with  all  other  atoms  with 
the  composite  potential  (Method  1);  or  they  could  be  treated  as 
foreign  defects  and  allowed  to  interact  with  only  a  repulsive 
potential  (Method  2).   A  total  of  three  different  defect  positions 
were  simulated:  an  interstitial  in  a  \010)  open  channel  (int  A), 
an  interstitial  in  a  < 100)  or  (001)  open  channel  (int  B) ,  and  a 
replacement  atom  (rep)  in  a  lattice  site.   (See  Figure  7.) 

1 .   Interst it ials 

As  previously  mentioned,  all  tungsten  interst it ials  initi- 
ated crovvdion  migration  when  treated  by  method  1,   An  interstitial 
treated  by  method  2  also  initiated  crowdion  migration  unless  buried 
in  the  center  of  an  enlarged  l4  x  l4  x  14  lattice.   Because  an 
interstitial  that  has  pushed  its  neighbors  away  has  a  lower  po- 
tential than  one  that  has  not  done  so,  the  numerical  values  for 
binding  energy  of  tungsten  interstitials  could  not  serve  as  a 
true  standard  for  comparison  with  W-Ne  results.   Qualitatively, 
however,  much  could  be  learned  from  the  W-W  energy  levels.   First, 
it  was  expected  that  all  energy  levels  of  a  defect  found  by 
Method  1  would  be  negative  at  equilibrium.   Values  for  the  compo- 
site potential  can  be  either  negative  or  positive,  but  are  posi- 
tive only  at  very  small  separations.   A  negative  potential  energy 
means  the  atom  is  bound  in  the  crystal.   As  shown  in  Figure  1,  on 
the  next  page,  the  first  two  interstitial  levels  for  the  W-W 
reaction,  Method  1,  were  at  -3-0  cV  and  -3.1  eV,  corresponding  to 
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the  two  inter st itials  in  the  first  layer  (int  A  and  int  B) .   The 
next  two  levels  were  at  -5.4  eV  and  -5.9  eV,  corresponding  to  the 
interstitials  in  the  second  layer.   Interstitials  in  the  third 
layer  and  deeper  had  binding  energies  ranging  from  -7.4  eV  to 
-8.1  eV.   The  value  of  -8.1  eV,  labeled  Mo°"  was  obtained  from  the 
interstitial  placed  in  the  center  of  the  seventh  layer  of  a 
14  x  14  x  l4  lattice.   Since  this  value,  and  all  other  values  for 
Method  1  binding  energies  were  reached  after  crowdion  migration, 
they  were  all  expected  to  be  lower  than  they  would  be  without 
crowdion  migration.   The  only  way  the  true  binding  energy  of  a 
tungsten  interstitial  could  have  been  found  would  have  been  to 
find  a  tungsten  crystal  size  large  enough  to  contain  the  crowdion 
migration  of  a  tungsten  interstitial  with  the  long-range  composite 
potential.   Compute  running  time  madp  this  impossible. 

Also  shown  in  Figure  1  are  the  binding  energies  for  the 
Method  2  W-W  interstitials.   Note,  first,  that  they  were  all 
positive.   This  was  again  expected,  since  the  potential  equation 
for  Method  2  is  positive  over  all  space.   Note,  second,  that  the 
binding  energies  for  interstitials  in  the  first  two  layers  were 
zero.   This  was  because  all  purely  repulsive  atoms  in  these  first 
layers  escaped  the  lattice  completely.   The  other  positive  levels 
shown,  were  incomplete,  because  many  interstitial  positions  did  not 
possess  stable  energy  levels.   The  unstable  levels  oscillated  be- 
cause of  significant  lattice  motion,  caused  by  crowdion  migration, 
and  measured  by  a  short  range  potential.   The  range  was  so  short, 
that  significant  jumps  in  binding  energy  occured  when  an  atom 
moved  into,  or  out  of  the  range  of  the  potential.   Evidently, 


33 


small  vibrations  in  atoms  with  an  equilibrium  distance  of  about 
^~3~  LU  from  the  interstitial,  frequently  caused  crossings  of  this 
range  limit,  adding  or  subtracting  energy  from  the  binding  energy 
each  time  one  of  them  crossed,  thus  invalidating  many  of  the  inter- 
stitial binding  energies. 

The  energy  levels  shown,  however,  demonstrated  the  meaning 
of  a  positive  "binding  energy".   The  numbers  reflected  the  "amount 
of  repulsion"  associated  with  different  positions  in  the  lattice. 
The  ordering  of  the  levels,  i.e.,  higher  energies  for  deeper  layers 
was  expected.   An  interstitial  deeper  in  the  lattice  felt  more 
repulsion,  because  it  was  surrounded  by  a  greater  number  of  re- 
pulsive neighbors.   Again,  the  level  labeled  "  °°"  represented  an 
interstitial  in  a  l4  x  l4  x  l4  lattice;  but  in  this  case  of  a 
Method  2  interstitial,  crnwdion  migration  did  not  occur. 

The  concept  of  a  positive  binding  energy  may  or  may  not 
be  the  actual  physical  situation,  but  it  is  still  academically 
valuable.   Instead  of  an  atom  resting  near  the  bottom  of  a  po- 
tential well,  as  in  Method  1,  an  atom  can  be  "wedged"  between  the 
repulsive  walls  of  its  neighbors.   In  both  cases,  the  atom  is 
"bound"  in  the  lattice.   The  ordering,  spacing,  and  other  cor- 
respondences between  the  positive  and  negative  levels  validate  the 
qualitative  use  of  the  positive  levels. 
2 .   Replacement  Impurities 

A  tungsten  lattice  with  a  replacement  impurity  is  a  per- 
fect tungsten  lattice,  but  Method  1  or  Method  2  could  be  used  on 
the  atom. in  question.   For  Method  1,  a  perfect  crystal  wasallowed 
to  relax  with  time,  yielding  only  negligible  motion;  and  the 
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binding  energies  of  the  center  atom  in  each  layer  was  recorded 
(see  Figure  1).   The  binding  energies  of  atoms  in  the  first, 
second,  and  third  layers  were-5.3  eV,-6.8  eV,  and-7.0  eV  respec- 
tively.  The  subsequent  reversal  of  order  for  levels  corresponding 
to  deeper  layers  is  a  program  anomally,  caused  by  the  use  of  a 
finite  depth  crystal.   Runs  on  larger  crystals  indicated  that  the 
order  would  not  reverse  in  an  infinite  lattice.   The  Mao"  level 
at~8.8  eV  is  the  experimentally  determined  heat  of  sublimation  [26] 
of  tungsten.   These  numerical   values  for  the  binding  energies 
were  valid  as  standards  of  comparison  for  the  Ne-W  data,  since  the 
motion  and  equilibrium  positions  of  the  replacement  atoms  and  their 
neighbors  were  nearly  identical,  and  usually  less  than  .3  LU  in 
both  cases.   Note  that  the  binding  energies  of  the  Method  1  re- 
placement atom??  wprp  lowpr  than  the  interstitial  atom  levels.   As 
previously  stated,  this  was  to  be  expected,  since  a  replacement 
atom  rests  in  the  bottom  of  a  periodic  potential  well  in  the 
lattice,  whereas  an  interstitial  rests  in  a  higher  well,  because 
it  is  nearer  to  its  neighbors  than  the  normal  equilibrium  sepa- 
ration.  Also  note  that  if  the  entire  Method  1  spectrum  of  binding 
energy  levels  were  used  as  a  standard  of  comparison  for  Ne-W  levels, 
then  the  W-W  interstitial  levels  should  be  higher  with  respect  to 
the  replacement  levels,  than  shown  in  Figure  1,  because  of  crowdion 
migration . 

The  Method  2  replacement  level  labeled  " a"'  cor  responded 
to  a  replacement  atom  in  the  center  of  the  fifth  layer  in  a  10  x 
10  x  10  lattice.   Note  that  it  was  not  above  the  interstitial 
levels,  as  it  should  have  been  by  comparison  with  Method  1.   This 
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was  the  major  weakness  of  Method  2:  because  it  was  a  measure  of 
repulsion  only,  and  neglected  the  potential  well,  it  under  esti- 
mated the  values  of  the  binding  energies  of  atoms  whose  normal 
position  was  in  that  well;  i.e.,  replacement  defects. 

E.   THE  NEON  DEFECT  IN  TUNGSTEN 

The  neon  atom  defect  was  again  placed  in  any  one  of  the  three 
lattice  positions,  labeled  "int  A",  "int  B" ,  and  "rep".   The  neon 
defect  could  be  treated  by  Method  2  only.   The  neon  energy  levels 
are  shown  in  Figure  2,  on  the  next  page.   Note  that  the  energies 
have  been  multiplied  by  a  mass  correction  factor.   (See  Appendix  D. ) 

1 .   Inter stitials 

The  neon  interstitial  never  initiated  crowdion  migration, 
but  as  in  the  W-W  case,  atoms  placed  in  the  first  two  layers  es- 
caped the  crystal,  and  therefore  had  zero  binding  energy.   Again, 
the  ordering  of  the  levels  was  a  measure  of  the  replusion  on  the 
interstitial,  which  increased  as  the  interstitial  was  placed 
deeper  in  the  lattice.   Again  the  M<:o"  level  was  the  result  of 
placing  a  neon  interstitial  in  the  center  of  a  l4  x  l4  x  l4  tung- 
sten lattice. 

2  .   Replacement  impurities 

Again  note  that  the  replacement  levels  are  lower  than 
would  be  expected  by  comparison  to  the  W-W  Method  1  standard.   It 
is  hypothisized  that  these  replacement  levels  should  be  higher 
than  the  interstitial  levels,  by  comparison  to  the  standard.   This 
assumption  is  valid  if  a  Ne-W  potential  well  exists.   It  is  fea- 
sible that  since  neon  is  almost  incapable  of  binding,  that  no  N'e-W 
well  exists.   On  the  other  hand,  the  ionization  state  of  both  noon 
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and  tungsten  in  the  lattice  is  unknown,  so  the  existence  of  a 
shallow  potential  well  is  quite  possible. 

F.   CORRELATION  WITH  EXPERIMENT 

1 .   Scaling  the  Levels  and  Peaks 

Since  no  direct  way  of  transforming  the  Ne-W  energy  levels 
into  correctly  scaled  negative  values  exists,  an  arbitrary  linear 
scaling  factor  between  KS's  data  and  the  Ne-W  levels  has  been 
used.   Note  that  every  level  or  group  of  levels  corresponded  to 
an  experimental  peak  in  Figure  2,  except  in  two  places:   first, 
the  broad  peak  at  about  2000   K  had  no  energy  level  counterpart, 
but  could  be  assumed  to  correspond  to  the  closely  ordered  replace- 
ment levels,  shifted  above  the  interstitial  levels  by  the  above 
hypothesis.   Second,  the  narrow  peak  at  450  K  was  without  an 
energy  level  counterpart.   Note  that  three  interstitial  levels 
had  zero  simulated  binding  energy  because  they  had  escaped  the 
crystal.   If,  however,  the  assumption  that  a  Ne-W  well  exists  was 
true,  then  some  or  all  of  these  three  interstitial  locations, 
(i.e.,  two  surface  layer  positions  and  the  open  channel  position 
in  the  second  layer)  would  be  stable,  bound  positions,  and  would 
be  expected  to  generate  an  energy  level  in  the  vacinity  of  the 
450   K  peak.   If  the  arbitrary  scaling  of  the  peaks  and  levels  was 
correct,  then  the  peak  at  450  K  is  proof  that  a  Ne-W  well  exists, 
since  a  purely  repulsive  potential  would  not  allow  a  first  layer 
or  second  layer  open  channel  interstitial  energy  level  to  exist. 
The  existence  of  this  well,  then,  would  in  turn  substantiate  the 
shift  of  the  replacement  levels  above  the  interstitial  levels,  to 
correspond  to  the  2000   K  peak. 
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2 .   Probes  of  Potential  Wells 

Various  attempts  to  substantiate  the  scaling  between  the 
levels  and  peaks  were  made.   No  approximation  to  a  Ne-W  potential 
well  could  be  justified,  and  the  correspondence  between  W-W  and 
Ne-W  results  was  not  complete  enough  to  invert  and  scale  the  po- 
tential energy  levels  to  realistic,  negative,  binding  energies. 

Attempts  were  also  made  to  match  both  Method  1  and  Method  2 
W-W  energy  levels  to  KS ' s  data.   The  Method  2  levels  were  too 
incomplete,  and  the  Method  1  levels  required  an  unknown  arbitrary 
reduction  of  interstitial  energy  levels  to  compensate  for  crowdion 
migration.   Too  many  alternate  reductions  were  possible  to  choose 
one  as  the  correct  rescaling. 

Another  approach  that  yielded  little  information  was  a 
plot  of  the  differences  between  interstitial  and  replacement 
energy  levels  for  each  layer  in  the  lattice. 

One  valuable  method  of  investigating  the  nature  of  possible 
Ne-W  negative  binding  energies  was  to  probe  the  perfect  lattice 
with  an  interstitial  at  various  initial  positions  and  plot  the  re- 
sultant potential  energy  of  an  interstitial  vs.  position.   Since 
the  potential  was  always  positive,  the  results  of  this  investigation 
were  potential  wells  above  the  x-axis .   Although  the  positive  lo- 
cation of  these  wells  was  not  realistic,  the  relative  depth  of  the 
wells  was  significant.   The  average  depth  of  a  neon  interstitial 
well,  deep  in  the  lattice,  was  about  4.2  eV  (see  Figure  8).   The 
graph  was  made  by  placing  inter stitials  in  positions  in  the  (010) 
open  channel,  and  thus  represents  the  barriers  that  the  interstitial 
must  penetrate  as  it  escapes  the  crystal.   Note  that  the  wells  v.o:e 
not  at  the  obvious  holes  in  the  BCC  lattice,  but  were  between  layers 
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In  the  actual  simulation,  the  interstitial  rarely  fell  into  this 
well,  but  instead  pushed  its  two  nearest  neighbors  away  and  made 
the  initial  position  the  low  potential  position.   Here  again,  sur- 
face effects  and  relaxation  reduced  the  tendency  for  inter st it ials 
to  relax  into  expected  infinite  lattice  equilibrium  positions. 
Slight  differences  in  the  final  equilibrium  positions  were  not  of 
significant  importance  to  binding  energies.   The  actual  numerical 
values  for  the  depths  of  these  positive  wells  were  not  necessarily 
scaled  properly  since  they  ignored  the  actual  Ne-W  potential  well, 
and  were  found  from  perfect  lattice  probes  at  time  zero,  before 
relaxation.   Nevertheless,  the  depth  of  a  well  deep  in  the  lattice 
of  4.2  eV  agreed  well  with  KS's  prediction  of  4-5  eV  [8]  for  the 
desorbtion  energy  corresponding  to  1720  K,  at  the  front  edge  of 
the  highest  pe^k.   1720  K    closely  approximated  the  temperature 
that  the  arbitrary  scaling  had  assigned  to  a  deep  interstitial. 

Note  that  the  initial  position  of  a  replacement  atom   was 
\[3~  LU  from  its  neighbors,  and  thus  because  of  potential  erosion, 
it  had  a  potential  of  zero  eV  initially.   To  climb  out  of  this 
well,  about  17-5  eV  must  be  supplied.   Although  this  number  was 
inaccurate  for  the  same  reasons  listed  above,  it  did  demonstrate 
that  even  a  purely  repulsive  potential  can  predict  a  greater 
binding  energy  for  replacement  atoms  than  for  interstitial  atoms. 
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V.   CONCLUSIONS 

An  arbitrary  scaling  has  been  used  to  correlate  the  simulation 
results  with  experimental  data.   Although  the  method  was  not  ana- 
lytically sound,  no  other  avenues  of  approach  to  the  problem  could 
be  found  that  could  further  justify  our  hypothesis. 

Satisfaction  can  be  gained,  however,  from  the  fact  that  these 
results  compare  favorably  with  known  data  at  many  interfaces.   Our 
model  was  a  tried  and  proven  one,  with  many  successful  sputtering, 
channelling,  and  similar  simulations  to  its  credit.   This  present 
model  invariably  behaved  in  a  physically  valid  manner  or  a  manner 
which  could  be  made  physically  acceptable  by  varying  the  con- 
trolling parameters  in  the  program.   Specifically,  many  previous 
experimental  and  simulated  results  for  infinite  crystals  were 
reproduced  when  simulation  took  place  deep  in  a  large  lattice, 
such  as  the  (lio)  split  interstitial  position  for  BCC  structures. 
The  Method  1  replacement  levels,  if  found  for  a  much  deeper  cry- 
stal, would  have  asymptotically  approached  very  close  to  the 
8,8  eV  heat  of  sublimation.   The  simulated  depth  of  the  positive 
potential  well  for  interst itials  of  about  4.2  eV  closely  approxi- 
mated KS's  prediction  of  4.5  eV. 

All  avenues  in  additional  computer  simulation  have  not  been 
exhausted.   Future  simulation  of  Argon,  Krypton,  and  Xenon  defects 
in  tungsten  should  be  fruitful.   Comparisons  between  the  relative 
locations  of  these  new  energy  levels  might  further  substantiate 
this  research.   In  particular,  if  simulation  can  explain  why  KS's 
neon  data  contains  five  desorbtion  peaks,  while  their  Argon, 
Krypton,  and  Xenon  data  contain  four  peaks,  it  will  be  a  major 
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success.  KS  also  gathered  data  for  different  crystal  surfaces 
and  different  angles  of  incidence,  which  might  be  investigated 
by  computer  simulation. 

This  simulation  was  also  important  in  that  it  investigated 
the  lattice  surface;  a  topic  which  has  rot  received  as  much 
attention  as  infinite  crystal  dynamics.   Since  radiation  damage 
theory  and  modern  transistor  theory  is  very  much  concerned  with 
the  crystal  surface,  the  computer  simulation  field  will  undoubtably 
increase  their  emphasis  on  surface  effects,  with  considerable 
attention  toward  better  ways  of  treating  a  foreign  interstitial. 

A  new  exhaustive  book  which  reports  on  the  present  state  of 
knowledge  in  all  these  areas,  with  emphasis  on  experimental  re- 
sults has  just  been  published.   It  is  a  report  on  the  proceedings 
of  the  International  Conference  on  Vacancies  and  Interstit ials  in 
Metals,  1968  [27'J  . 


42 


APPENDIX  A:  CRYSTAL  GEOMETRY 

The  computer  program  can  call  any  one  of  nine  lattice  generator 
subroutines:  three  face-centered  cubic  subroutines  (LlIOO] ,  [lIIO], 
[LIU]),  three  body-centered  cubic  subroutines  ([bIOO],  [bIIO], 
[Bill]),  and  three  diamond  subroutines  ([dIOO],  CdIIO],  [Dill]). 
The  diamond  subroutines  are  never  used  and  therefore  not  compiled; 
but  provision  has  been  made  for  their  future  inclusion  in  the  pro- 
gram.  The  dimensions  of  the  lattice  chosen  were  controlled  by  the 
input  data  variables  [ix],  Liy],  and  [ iz] .   Each  atom  in  the  cry- 
stal was  numbered,  in  the  order  x  followed  by  z,  followed  by  y. 
For  the  surface  layer  (Y  =  0)  of  the  tungsten  10  x  10  x  10  lattice, 
atoms  were  numbered  from  2-26;  for  the  first  layer  below  the  sur- 
face, atoms  were  numbered  27-51,  etc.   Atom  number  1  was  the  pri- 
mary, 01  point  defect  atom.   lLDj  was  the  ijumber  of  the  last  mobile 
atom,  or  the  last  atom  in  the  eighth  layer,  number  200;  and  LllJ 
was  the  number  of  the  last  atom  in  the  crystal,  number  250. 

The  placement  of  point  defects  was  accomplished  as  follows: 
after  the  desired  perfect  lattice  was  built  to  the  desired  size, 
subroutine  PLACE  was  called.   Three  types  of  defects  were  allowed: 
vacancies,  inter st it ials ,  and  replacement  impurities.   The  type 
and  location  of  the  defect  were  controlled  by  input  data  vari- 
ables:  the  type  of  defect  [iTYPe],  an  atom  number  [nVAC]  and  a 
displacement  vector  CdIX,  D1Y,  Dlz]  in  LU .   If  ClTYPEj  =  1,  a 
vacancy  was  created  in  site  number  [NVAC] .   This  "removal"  was 
accomplished  by  setting  ClCUT  (NVAC)]  =  1  which  "turned  off"  the 
atom,  removing  it  from  all  calculations.   If  LiTYPE]  =  2,  an 
interstitial  was  created  in  a  position  [-DIX,  -DIY,  +  Dlz]  LU  from 
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site  number  lNVACj .   This  interstitial  was  always  atom  number  1. 
Since  atom  number  2  was  always  at  the  origin,  number  1  could  be 
placed  using  a  displacement  vector  from  the  origin  (from  [nvac]  =  2) 
or  using  a  displacement  vector  from  a  site  next  to  the  interstitial. 
If  L ITYPe]  =  3,  a  vacancy  was  created  in  site  number  LNVAC]  and 
a  replacement  impurity,  put  in  its  place.   Note  that  for  both 
[iTYPEj  =  2  and  3,  either  a  foreign  or  self  defect  could  be 
placed.   For  the  case  of  either  the  self -inter stitial  or  the  self- 
replacement  atom  (giving  us  back  the  perfect  crystal),  either 
method  1  or  method  2  of  calculating  the  potential  could  be  used. 
The  choice  of  methods  was  also  an  input  parameter:  for  method  1, 
[iq]  =  1,  and  for  method  2,  [ IQ]  =  2. 
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APPENDIX    B:    POTENTIALS 

(In    this    appendix   and   all    subsequent    appendices,    the   brackets 
denoting   program    language  are   dropped.      Program    language    is    still 
written    in   all    capital    letters.) 
A.  BORN-MAYER    REPULSIVE   POTENTIAL 

1.  Potential    Energy:    For    the    lattice   atom    interactions,    the 
Born  Mayer    potential    equation    is: 

V.  .    =    exp(A    +    Br  .  .)  (1) 

°r  POT    =    EXP(EXA    +    EXB*DIST)  (1A) 

For    bullet-lattice   atom    interactions, 

V.  .    =    exp(A'    +    B'r  .  .)    -   V.  .(ROE)  (3) 

where   V.  .(ROE)    is    subtracted   to   retard   the  potential    so   that    it 

goe^    to   zero   at    the   nearest    neighbor    distance.       in    the  progiam, 

the  V.  .    equation    is: 
ij 

POT    =    EXP(PEXA    +    PEXB*DIST)     -    PPTC  (3A) 

2.  Force:    For    the    lattice  atom    interaction, 

-bv.  . 

Force   =  t — ^  =    -B   exp(A   +   Br  .  . ) 
Or  .  .  r v  ij 

=   exp    [(^2   -B   +   A)    +    Br  .  .]  (4) 

in    the   program,   0n(-B   +   A)     =    ( ALOC-C -EXB*CVEd] +    EXA)    =    FXA, 
where  CVED   is    a    conversion    factor    for    units,    and 

FORCE    =    [FXA    +    EXB*DIST]  (4A) 

For    bullet-lattice  atom   forces, 

-0V.  . 

Force    =  -s—^1  =    -    B' exp(  A' +B '  r  .  .)    =    expOi-B'+A1  )+3'r  .  .]    (5) 

Or  .  .  r v  in'  ij 

where    (^-B1  +A*  )  =  (  ALOC-C  -PEXB*CVED] +PEXA )    =PFXA,    and 

FORCE    =    EXP[PFXA    +    PEXB*DIST]  (5A) 
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=    exp[0n(2a)+(0nD+2ar    )  -  (  2a  )r  .  .]-exp!>ia  +(&n(  2D)  +ar    )-(a)r..] 

=  exp[alog( 2*alpha*cvr*cved) +alog( dcon) +2 . *alpha*re 

- ( 2 . *alpha*cvr ) *dist] 

-  exp[ alog( alpha*cvr*cved) +alog( 2 . *bcon) +alpha*re 

-( alpiia*cvr)  *dist'j 

=  exp[ alog( -cgb1*cved) +cgd1 ) -cgb1*dist] 

-  exp[(alog(-cgb2*cved)+cgd2)-cgb2*dist"j 

=  exp[cgfi-cgbi*dist]-  exp[cgf2-cgb2*dist]  (7a) 

c.     cubic  fit 

1.  Potential  Energy:  The  best  cubic  fit  between  the  BM  and 
Morse  potentials  is  calculated  in  Subroutine   CROSYM.   The  po- 
tential equation,  defined  between  ROEA  and  ROEB,  is 

?  2 

POT  =  CP3r  .  .   +  CP2r  .  .   +  CPlr  .  .  +  CP0       (S) 
ij         iJ         iJ 

or  POT  +  DIST*(DIST*(DIST*CP3+CP2)+CP1)+CP0     (8a) 

-&POT  2 

2.  Force:  Force  =   x     =  -3CP3* •  .  -2CP2r .  .  -  CP1    (9) 

Or  .  .  in        ii 

lj 

2 
=  (-3-*CP3*CVED)r . .  +  (-2.*CP2*CVED)r  .  .+  (-CPl*CVED) 

2 
=  CF2r  .  .   +  CFlr  .  .  +  CF0 ,  or 

FORCE  =  DIST*(DIST*CF2+CFl)+CF0.  (9A) 
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APPENDIX  C:   AVERAGE  FORCE  METHOD  AND  TIME  DURATION  THEORY 

A.   AVERAGE  FORCE  METHOD:  The  average  force  technique  has  been 
explained  in  great  detail  in  Ref.  15.   It  was  summarized  in 
Chapter  3,  and  therefore  discussion  here  is  limited  to  the  average 
force  method  in  the  program  language. 

When  the  desired  lattice  is  built,  the  position  of  the  ith  atom  is 
stored  simultaneously  in  RX(I),  RY(I),  RZ(I);  RXK(I),  RYK(I), 
RZK(I);  and  RXI(I),  RYI(I),  RZI( I) .   The  latter  set  of  coordinates 
never  change  and  are  used  for  comparing  new  positions  to  original 
positions,  and  in  calculating  DX(I),  DY(I)  and  DZ(I)  for  output. 
The  middle  set  of  coordinates  containing  the  letter  K  are  for 
storing  the  initial  positions  at  the  beginning  of  each  timestep. 
A  step  by  step  summary  of  the  average  force  method,  showing  X 
coordinate  calculations  only,  follows: 

1.  Based  on  the  position,  RX(I),  of  the  ith  particle  at  the  be- 
ginning of  the  timestep,  the  force  FX(I)  is  calculated  in  STEP. 

2.  RX(I)  is  stored  in  RXK(I). 

3.  The  new,  temporary  position,  RX(I)  is  calculated,  based  on 
the  force  at  RXK(I) : 

X1  =  X  +  VAt  +  F(At)  /2m  (10) 

or 

RX(I)  =  RX( I)+DTOD*(HDTOM*FX(I)+VX( I) .     (10A) 

4.  A  new  force  is  calculated,  based  on  this  new  position  RX(I). 

5.  VSS  stores  VX(I),  the  original  velocity  of  the  ith  atom.   This 
velocity  is  half  the  velocity  of  the  ith  atom  in  the  previous 
timestep :  the  \   factor  being  an  arbitrary  damping  multiplier.   A 
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new  velocity,  based  on  the  new  force,  is  found: 

V1  =  V  +  FAt/2ra  (11) 

or 

VX(I)    =    VSS+HDTOM*FX(  I)  .  (HA) 

6.      The   final    position    is    calculated,    based   on   the   average   of 

these   two   velocities 

X1    =   X   +   %At(V+V1)  (12) 

RX(I)     =    RXK(I)+(VX(I)+VSS)*HDTOD  (12A) 

The  resultant    velocities    are   halved,    a    new    timestep    duration    is 
calculated,    and   the   process    repeated. 

B.       TIMESTEP    DURATION    THEORY 

This    simulation    uses    the   best   possible   estimate    of  a    timestep 
duration,    DT,    as    calculated   from   the   present    state    of   the    forces 
and    energies     in    the    lattice   for    use    in    the    next    timestep.       To 
limit   motion    to  an    increment    small    enough    to  preserve   the   accuracy 
of    the  average    force   appr  oximation ,    we   define    DTI   as    the   maximum 
distance  any    atom    is   allowed   to  move    in    one   timestep. 
From    (10),    we    find 


Therefore, 


For 


AX.    =    (V.+    F.At/2m)At. 
l  v     l         l  ' 


At    =    AX./(V.+F.At/2m) .  (13) 


V.    »    F.At/2m,       At    =    AX./V..  (14) 

li  ll 


If   we   find    the   fastest    moving   atom   and  assure    that    it    does    not 
move  more    than    DTI,    we    have    limited   the   motion    of   all    other   atoms 
to    less    than    DTI. 

ThUS'  DT    =    (DTI*CVD)/EMAX    =    FDTI/EMAX  (l^A) 
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where 

EMAX    =    SQRT(VX(I)*VS(I)+VY(I)*VY(I)+VZ( I)*VZ( I) . 

For  V.    »    F.At/2m, 

1  1 

2raAX. 
At    =    AX./F.At/2m    =      — —       .  (15) 

l 

Anatagous  to  above,  we  find  the  most  stressed  atom  and  assure 
that  it  does  not  move  more  than  DTI.   Thus, 

DT  =  SQRTC  (2.*PTMAS*DTI*CVI))/FMAX] 

=  sqrt[tfac/fmax]  (15A) 

where      FMAX=    SQRTC  FX(  I )  *FX(  I)  +FY(  I)  *FY(  I)  +  FZ(  I)  *FZ(  I)  "J  . 
Since   rigorously   we   cannot    make   either    or    these    limiting   assumpt- 
ions,   we  must    go   back   to   our    original    equation    for    DT ,    equation    (13) 
Since   this    equation    involves    DT ,    we  proceed  as    follows; 

1.  Assume   V.    «    F.At/'2m   and   calculate    At    from    (15a). 

2.  Insert    this   preliminary    value   for    DT    in    (13)    and   compare 

V.    to    F.At/2m.       If    V.     is    larger,    calculate   DT    from    (lAA).       If 
11  l  v  ' 

F.At/2m    is    larger    calculate    DT  from    (15A). 

A   complication   arises    when    a  foreign    impurity    is    in    the    lattice, 
because   of  a    variation    in    the   value    for    m    in    (15).      This    is 
especially   acute   when    the    differences    in   masses  are   great.      The 
method   used   to   solve   this   problem    is    as    follows: 

1.  If    either    FMAx   =   F      or_  EM  A  X  =   V    ,    the   entire   proceedure,    above, 
is    followed  using    the   mass    of   the  bullet    for    m. 

2.  If   both   FMAX   /   F_    and   EMAX   f  V    ,    the    entire  proceedure    is    fol- 
lowed  using    the   mass    of   a    lattice  atom. 

The    requirement    that    the   bullet   mass    be    used   if   either    EMAX   or 
FMAX  describe   the   bullet    circumvents    the  problem    of    having    the 
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bullet  the  fastest  moving  atom,  but  not  the  most  stressed  atom, 

or  visa  versa. 

-14 
To  begin  the  problem,  an  arbitrary  value  of  10     seconds  is 

assigned  to  DT.   If  at  any  time  in  the  program  EMAX  =  FMAX  =  zero, 

-14  ...  ... 

10     is  again  assigned  to  DT  to  prevent  division  by  zero. 
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APPENDIX  D:  SUBROUTINES  STEP,  ENERGY,  AND  LOCAL 

A.  DISTANCE  CALCULATIONS 

In  all  three  subroutines,  STEP,  ENERGY,  and  LOCAL,  a  method  of 
finding  all  atoms  within  a  given  radius  of  another  atom  was  needed. 
For  lattice  atom  interactions,  atoms  inside  ROEA,  ROEB,  and  ROEC 
were  found;  for  the  foreign  interstitial  interactions,  atoms  inside 
ROE  were  found;  and  for  LOCAL,  atoms  inside  ROEL  of  a  point  defect 
were  found.   The  time  saving  technique  used  to  do  this  was  to 
successively  eliminate  all  atoms  with  an  x  component  difference 
greater  than  the  given  radius,  then  similarly  for  y  components, 
then  for  z  components.   The  resulting  volume  not  eliminated  is  a 
cube  circumscribing  the  desired  sphere.   Finally  the  time-con- 
suming test  of  eliminating  all  atoms  for  which  the  desired  radius 
is  less  than  SQRT  ( ERX*DRX  +  DRY*DRY  h  DRZ*DRZ)  is  applied  to 
only  atoms  inside  the  cube. 

B.  SUMMATION  INDICES  IP  AND  IQ 

Interactions  V.  .,  0.  .,  and  F.  .  are  found  by  evaluating  all  values 

in  the  half  matrix.   For  example  F   ,  F   ,  F  ,  •••  are  found,  then 

F„_,  F_i  ,  •••  etc.   The  variable  IP  controls  the  starting  point 
23    3h 

for  the  j  summation.   IP  is  always  set  to  I  +  1  to  avoid  the  re- 
petition of  finding  F..  and  F...   IQ  controls  the  starting  point 
for  the  i  summation.   If  the  primary  is  to  be  treated  as  a  lattice 
atom,  IQ  =  1.   If  it  is  to  be  treated  as  a  foreign  particle, 
IQ  =  2,  and  all  F.  .  are  found  separately. 

C.  DISTRIBUTION  OF  FORCES  AND  POTENTIAL  ENERGIES 

The   forces    F..   are   equal   and    opposite   on    i   and    j:       i.e.,    F.=-F.. 
The  potential    energies   are   split    in   proportion    to   the   reduced 
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mass  of  the  interacting  particles.   This  is  easily  understood  by 
observing  the  kinetic  energy  distribution  of  an  elastic  collision 
of  m  and  M  where  M  >   m.   We  find  that  m  carries  away  almost  all 
the  kinetic  energy:  specifically  it  carries  away  I — — 1  E   ~E    . 
If  a  pair  of  atoms  are  to  behave  elastically,  the  potential 
energies  which  are  transformed  into  kinetic  energies  of  motion 
must  be  split  in  the  same  manner.   For  this  reason, 

PPE<1»=  (toas^as)  *  POT  ■  BSAVE*POT; 

and 

PPE<J>=  (tMAS^Ias)    *P0T    =   TSAVE*POT. 

note,    for    BMAS   =   TMAS,    BSAVE   =   TSAVE   =  h,    and   the   energies    are 
split    equally. 

D.       SUBROUTINE   LOCAL 

LOCAL  measures    the   change    in  potential    energy   associated   with  a 
sphere   of   radius   ROEL   surrounding    a   point    defect.       It    sums    up   the 
potential    energies    of    each  atom   found    inside    this  sphere   at    time 
zero.       It    remembers    these   atoms,    and    for    each    timestep    re-sums 
the  potential    energies    of    these    same   atoms.      The   sum,    total    local 
potential    energy  TLPE,    is    subtracted   from  TLPE   at    time    zero 
(TLPE0) ,    to   give   a   measure    of   the   change    in   potential    energy 
(DLPE)    inside   ROEL. 
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APPENDIX   E:     COMPUTER   PROGRAM    GLOSSARY 

NOTE:       In    this    glossary,    the    terms    "point    defect    atom",    "bullet", 

and    "primary"   are   synonymous;    and   the   terms    "lattice  atom"   and 

"target"   are   synonymous. 

ALPHA:       Input    Morse  potential   parameter 

BSAVE:       Target    mass/(target    mass    +   bullet   mass);    distributes 

potential    energy   between    target    and  bullet 
BIND:         Negative    of   the   total   potential    energy    (TPOT)    at    time   zero 
BMAS:         Mass    of    bullet    in   amu 

BULLET:    Alpha-numeric   array    for   point    defect   material 
CFO,    CF1 ,    CF2:       Force  parameters    of   cubic    fit   between  Morse   and 

Born-Mayer    functions 
CGBl ,    CGB2:      Morse  potential   parameters 
CGD1 ,    CGU2:      Morse  potential   parameters 
CGF1 ,    CGF2:       Morse   force  parameters 
CPO,    CP1,    CP2,    CP3:    Potential   parameters    of    cubic    fit   between 

Morse  and   Born-Mayer    functions 

CVD:    CVR    x   10         ,    converts    lattice   units    to  meters 

-19 
CVE:    1.6   x    10         ,    converts    electron    volts    to    joules 

CVED:    CVE/CVD,    a    ratio   used   to  avoid   repeated   division 

-27 
CVM:    1.672    x   10         ,    converts    atomic   mass    units    to  kilograms 

CVR:    LU    in    angstroms;    converts    lattice  units    to  angstrom    units 

D1X,    D1Y,    DlZ:    Displacement    coordinates    for    location    of    interstitial 

from   reference  atom  ,      NVAC 
DCON:         Input   Morse  potential   parameter 
DFF:  ROE-DIST,    the   distance   closer    than   ROE    that    an    atom    is    to 

the   primary 


54 


DIST:         Distance   between   any    two   atoms 

DLPE:         TLPE-TLPE0,    the    change    in    total    local    potential    energy 

since   time    zero 
DRX,    DRY,     DRZ:    x,y,z    components    of    DIST 
DT:  Length    of   a    tiraestep    in    seconds 

DTI:  Number    of    lattice   units    most    energetic   atom   may   move    in 

one   times tep 
DTOD:         DT/CVD — a    ratio   used   to   avoid    repeated   division 
DTOM:         DT/PTMAS--a    ratio   used   to  avoid   repeated    division 
DTOMB:       DT/PEMAS--a    ratio   used  to  avoid   repeated    division 
DX(I),     DY(I),    DZ(I):    Change    in   position    of    it_h    atom    from    initial 

position   at    time    zero 
EMAX:         The   maximum    energy    encountered    in   any    cycle 
EV:  Primary    energy    in    electron    volts 

EVR :  Primary    energy    in    kilo-electron    volts 

EXA ,    EXB:    Input    Born-Mayer    potential    function   parameters    for    the 

target 
F2 :  Square   of    the  force    on   a    specific   atom 

FA:  The    component    force    increment    on   an    atom 

FDTI :    DTI    X   CVD,    a   parameter    used    to   determine    DT  my   maximum 

energy   method 
FM:  A   small    number    used   in   checking  potential    energy    zero 

point 
FM2:  FM   squared 

FMAX:         Maximum    total    force    on    the   most    stressed   atom    in    the 

crystal 
FOD:  FORCE/DIST--a    ratio   used   to   avoid   repeated    division 
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FORCE: 

FX(I), 

FXA: 

HBMAS : 

HDTOD: 

HDTOM: 

HDTOMB 

HTMAS: 

II: 

13: 

I DEEP: 

IH1 

IH2 

IHB: 

IHS: 

IHT: 

I  LAY: 

IN: 

IP: 

IQ: 


Numerical    value    of   the   force   function   with    a    variable 
parameter 
FY(I),    FZ(I):    x,y,z    components    of    total    force   on   an   atom 
Born-Mayer    force   function  parameter 
\   BMAS--a    ratio   used   to  avoid   repeated   division 
^§   DTOD--a    ratio   used   to   avoid  repeated   division 
\   DTOM--a    ratio   used   to  avoid   repeated   division 
%   DTOMB—a   ratio   used  to  avoid   repeated   division 
^  TMAS--a    ratio   used   to  avoid   repeated   division 
Variable    in    cubic    fit    subroutine 


11         it 


I SHUT: 

IT: 

ITT: 


Number    of    mobile    layers 
Alpha    numeric   array    for    program    title 
"  "  "  "        Morse   function   parameters 

"  "  "  "        bullet    element 

"  "  "  "         type  and    orientation    of    crystal 

"  "  "  "         target    element 

Same    as    I DEEP 

Odd-even    integer    used   to    determine   atom    site   establishment 
Subscript    value    of  atom.      Used    in   subroutines    STEP    and 
ENERGY 

Parameter    that    determines    whether    or    not   a    self   defect    is 
to   be  given   a    repulsive  potential    or    a    composite  attractive- 
repulsive  potential 

A  parameter    used   to   shut    down    the  program 

Unsealed   fixed  point    x    coordinate   used    in    lattice  generation 
Odd-even    integer    used   to   determine  atom   site    establishment 


56 


ITYPE:   Parameter  used  to  determine  the  type  of  point  defect: 

vacancy,  interstitial,  or  replacement 
IX,  IY,  IZ:  Number  of  x,y,z  planes  of  crystal 
J2 :      Variable  in  the  cubic  fit  subroutine 

JJ:      Parameter  in  the  BCC(lll)  lattice  generation  subroutine 
JT:      Unsealed  y  coordinate  used  in  crystal  generation 
JTS:     Variable  used  to  establish  atom  sites 
JTT :     "         "     "   "  "     " 

KF:      Final  K  in  LOCAT  (K)  assigned  to  an  atom 
KT:      Unsealed  z  coordinate  used  to  establish  atom  site 
LCUT( I) :  Used  to  identify  an  ith  atom  which  is  not  included  in 

calculations 
LD:      The  highest  numbered  atom  in  the  mobile  layers 
LL:      The  highest  numbered  atom  in  the  entire  crystal 
LOCAT(K) :  Dimensioned  variable  that  remembers  the  numbers  of  the 

atoms  within  a  radius  ROEL  of  the  primary  at  time  zero 
LS:      Variable  associated  with  each  of  the  nine  lattice 

generator  subroutines 
MCRO:    One  number  higher  than  the  order  of  the  fit  between  the 

Born-Mayer  and  Morse  potentials,  always  4  in  this  simu- 
lation 
ND:      Data  output  increment,  in  numbers  of  timesteps 
NEW:     Parameter  used  to  determine  whether  or  not  atom  numbers 

have  been  stored  in  LOCAT(K) 
NPAGE:   Page  numbering  variable 
NRUN:    Parameter  used  to  determine  whether  or  not  to  read 

additional  data  cards 
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NS:  Initial   print    statement    timestep    number 

NT:  Timestep    number 

NTT:  Timestep    number    limit    before   shutdown 

NVAC:         An   atom   number    used   to   establish   point    defects    or    used  as 

a    reference  point    for    interstitial  placement 
PAC:  Parameter    for    bullet    force   function    correction 

PBMAS:       Primary   mass    in   kilograms 
PEXA,    PEXB:    Input    Born-Mayer    potential    function  parameters    for 

the  bullet-target    interaction 
PFPTC:      Primary    force   function    evaluated  at    ROE 
PFXA:         Primary    force   function   parameter 
PKE(I):    Kinetic    energy    of    the    ith    atom 
PLANE:       Alpha-numeric   array    for    lattice   orientation 
POT: 

PPE(I):    Potential    energy    of   the    ith   atom 
PPTC:         Primary   potential    function    evaluated  at    ROE 
PTE(I):    Total    energy    of    the    ith    atom    (potential    +   kinetic) 
PTMAS:       Target    mass    in    kilograms 
RE:  Input   Morse  potential   parameter 

RO:  Spacing    constant    in    FCC(llO)    lattice   generation    subroutine 

ROE:  Nearest    neighbor    distance 

R0E2:         ROE   squared 

ROEA:  Maximum    cut    off    for    Born-Mayer    potential 

ROEB:         Minimum    cut    off   for    Morse  potential 
ROEC:         Maximum   cut    off   for    Morse  potential 
R0EC2:       ROEC    squared 
ROEL:         Radius    inside   of   which    local  potential    energy    is    found 


Potential    energy   between    two   atoms 
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ROEL2:       ROEL    squared 

ROEM:         ROE-DTI,    region    in   which   modification    of   repulsive   force 

must    be  made 
RX(I),    RY(I),    RZ(I):    x,y,z    coordinates    of   an    ith   atom  at   any    time 
RXI(I),    RYI(I),    RZI(I):    x,y,z    coordinates    of   an    ith   atom's    initial 

position 
RXK(I),    RYK(I),    RZK(I):    x,y,z    coordinates    of    temporary   position    of 

an    ith   atom    during    force   cycle 
SAVE:         h  POT 

SCX,    SCY,    SCZ:    x,y,z    coordinate   scale   factors 
SSCZ:         A   z    scale   factor    used   for    the   FCC(lll)    lattice   generator 

subroutine 
START:       An    optional    timing   variable,    not    used    in    this    simulation 
SUM:  Variable    in    cubic    fit    subroutine 

TARGET:  Alpha-numeric   array    for    target   material 

TSAVE:       Bullet   mass/( target   mass    +   bullet    mass);    distributes 

potential    energy   between   target    and  bullet 
TE:  Total    energy    of   all   crystal   atoms    (kinetic    +   potential) 

TEMP:         Temperature   of    lattice    in    degrees   Kelvin.      Not    used    in 

this    simulation 
TFAC:         A   time    factor    ratio   used   to   determine    DT   by  maximum   force 

method 
TFACB:       TFAC    for    the   bullet 

THERM:       Thermal    energy    of   atom.       Not    used    in    this    simulation 
TIME:         Elapsed  problem   time    in   seconds 

TLPE:         Total    local   potential    energy    of   atoms    within   a    radius    ROEL 
TLPE0:       TLPE   at    time    zero 
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TMAS:    Target  atom  mass  in  amu 

TPKE:    Total  kinetic  energy  of  all  crystal  atoms 

TPOT:    Total  potential  energy  of  all  crystal  atoms 

VSS:     Storage  variable  for  velocity  components 

VX(I),  VY(I),  VZ(I):  x,y,z  components  of  ith  atoms  velocity 

X,  Y,  Z:  Unsealed  coordinates  used  in  crystal  generation 

YLAX(I):  Relaxation  in  -y  direction  of  it_h  layer  in  L.U. 

ZP:      Floating  point  form  of  JTT 
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FIGUPE    7- 

BCC     (100)      ORIENTATION 
HARD     SPHERE     MODEL 
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FIGURE    9. 
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ERODED     POTENTIAL 
=  V*(R) 


ROEA 
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c 

THIS  PROGRAM  GENERATES  VARIOUS  TYPES  AND  ORIENTATIONS  OF 
CRYSTAL  LATTICESt  AND  INJECTS  A  VACANCY,  INTERSTITIAL,  OR 
REPLACEMENT  IMPURITY  AT  A  DESIRED  LOCATION,   IT  THEN,  BY  USE 
OF  ATOMIC  POTENTIAL  PARAMETERS  AND  NEWTONIAN  MECHANICS,  CAL- 
CULATES THE  DYNAMIC  RESULTS  OF  THE  SYSTEM;  OUTPUTING  POS- 
ITION, VELOCITY,  AND  ENERGY  VALUES  FOR  EACH  ATOM  IN  THE 
CRYSTAL. 
C 
DIMENSIONING  OF  VARIABLES  NOT  NEEDED  IN  COMMON 

DIMENSION  VX(IOCO)  ,  VY(  100C)  ,VZ(  10()o  )  ,  PKE(  10CC) 

DIMENSION  DX(IOOC) , DY( 1000) , DZ( 1000 ) , PTE( 1000) 

DIMENSION  RXK( 1000 ) ,RYK(1000 ) ,RZK( 10C0) 
C 
CCMMON  LABELING  OF  VARIABLES  REQUIRED  IN  OTHER  SUBROUTINES 

C0MM0N/CCM1/RX ( ltOC ) ,RY(1000 ) ,RZ( 100C ) »  LCUT( 1000 ) , 
1LL,LD, ITYPF,NVAC 

COMMON/COM2/IHK20)  , I H2 ( 8 ) , I HS ( 10 ) , IHB ( 6 ) , I HT ( 6 ) , 
1 TARGET (4) ,TMAS , BULLET (4 ) , BMAS , PLANE , TEMP , T HERM 

C0MM0N/C0M3/RXI ( lOoC) ,RYI( 1000) ,RZI( 1000) ,CVR,EVR, 
INT, TIME, DT, DTI , ILAY 

COMMON /C0M4/IX, I Y, I Z , SCX , SCY , SC Z , IDEE P , Dl X , Dl Y , Dl Z 

COMMON/COM5/ROE,ROE2,ROEM,EXA,EXB,PEXA,PEXB,FXA,PFXA, 
1 IQ,TSAVE,BSAVE 

COMMON/COM6/FX(1000) , F Y ( 10  00 ) , F Z ( 1000) ,PAC,PFPTC, FM 

COMMON /C 0M7/ PPT C,T POT, PPE( 1000 ) , TLPE , ROEL , R0EL2  ,  NEW 

C0MM0N/C0M8/R0EA,R0EB, ROEC , R0EC2 , CPO , CPl ,CP2,CP3, 
1CF0,CF1,CF2,CGD1,CGD2,CGB1,CGB2,CGF1,CGF2 

COMMON/CCMA/  A(4,5),MCR0 
C 

READ  STATEMENT  FORMATS 
9010  FORMAT! 20 A4) 

FORMAT (8A4,3F8. 5, 2F5.2) 

F0RMAT(4A4,3F8.5,6A4,F6.2) 

FORMAT (F6. 2,F5o3,I5,6I4,3F5.2, 12) 

FORMAT!  1UA4.A4.4I3  ,F8.4,I4) 


90  20 

90  30 

90  40 

90  50 

r 

o 

WRITE 

9610 

9620 

STATEMENT  FORMATS 

FORMAT! 1H1 ) 

F0RMAK47X, 'SUMMARY  OF  ATOMS  «//,  35X  ,  8A4  ,  •  ,  NT  =«I4,//, 
13!  '  ATOM       POSITION       BIND  ENERGY    '),//) 
96  30  FORMAT (3(  I  5 , 3F6. 2 , F 8. 4 , 8X )  ) 

9640  F0RMAT(/4X,F1C.3,25H  EV, TOTAL  KINETIC  ENERGY, , F10. 3 , 
127H  EV, TOTAL  POTENTIAL  ENERGY , FIG . 3 , 13H  EV , REDUCT I  ON, 
1//,20X,F10.3,50HEV,  LOCAL  POTENTIAL  ENERGY,  IN  VOLUME 
10F  RADIUS  =  ,F5.2,  /30X,  16HCHANGE  IN  TLPE  =  ,  FIG. 3) 
96  50  FORMAT! 10  5X,4HPAGE , 13, /,1H1) 
9660  FORMAT!/   ■       ATOM     DX         DY         DZ 

1VX         VY         VZ        KE         PE         TE«/) 
9670  FORMAT! 118, 3F10. 3, 3F10. 1,3F10,4  ) 
C 
INITIALIZING 

START=G.01*ITIME(XX) 

DO  2  1=1,1000 

LCUT! I )=C 

RXK! I }=0.0 

RYK! I )=0.0 

RZK! I )=OeU 

VX( I )=G.0 

VY(  I) =0.0 

VZ! I) =0.0 

PKE! I )=0.0 

PPE! I)=0.0 

PTE(I)=0.0 

RX( I »=C.O 

RY(I)=C.O 

RZ( I) =0.0 

FX( I )=0.0 

FY! I )=0.C 

FZ( I ) =0.0 

RXI ( I )=O.0 

RYK  I)=0.0 
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c 

INPUT 


RZI( I )=0.0 

ISHUT=1 

NRUN=0 


C 

CONSTA 


DATA 

READ 

READ 

READ 

READ 

READ 
DlXtDlY 

READ    ( 


5,9C1P  ) 
5,9C20  ) 
5,9030  ) 
5,9C30) 
5,9C40) 
D1Z,  IQ 
5,9050) 


IH1 

IH2,DCON,ALPHA,RE,ROEC,ROEL 

BULLET,BMAS,PEXA,PEXB,IHB,THERM 

TARGET, TM AS f EXA,EXB,1HT,TEMP 

EVR, DT I ,NTT,NS, ND, IP, IDEE P , I  TYPE , NVAC , 

IHS,PLANE,LS, IX, I Y, IZ,CVR,MCRO 


C 

REPULS 


NTS  AND  SCALING  FACTORS 

ROE2=3.0 

ROE=SQRT(ROE21 

ROEM    =    RGE-DTI 

R0EL2=R0EL*R0EL 

CVE=1.60E-19 

EV=EVR*1.0E+3 

CVM=1.672E-27 

FM=l„OE-10 

FM2=FM*FM 

CVD=CVR*l«OE-10 

CVED=CVE/CVD 

PTMAS=TMAS*CVM 

PBMAS=BMAS*CVM 

HTMAS=0.5*PTMAS/CVE 

HBMAS=U.5*PBMAS/CVE 

TSAVE=BMAS/( BMAS+TMAS) 

BSAVE=TMAS/ ( BMAS+TMAS ) 


C 

ATTRAC 


C 

CUTOFF 


C 
PARAME 

BETWEE 
(ROEA) 
TIAL  ( 
F  I T  T I  N 


IVE    POTENTIAL     PARAMETERS 
FXA=ALOG(-EXB*CVED)+EXA 
PFXA=ALOG(-PEXB*CVED)+PEXA 
PPTC=EXP(PEXA+PEXB*ROE) 
PAC~ ALQG(C^pn>  +dp y  a 
PFPTC=EXP( PAC+PEXB*ROE) 

TIVE    POTENTIAL     PARAMETERS 

CGDl=ALOG(DCON)+2.0*ALPHA*RE 

CGD2=ALOG(2.0*DCON)+ALPHA*RE 

CGB1=-2.0*ALPHA*CVR 

CGB2=-ALPHA*CVR 

CGF1=AL0G(-CGB1*CVED)+CGD1 

CGF2=ALOG(-CGB2*CVED)+CGD2 

DISTANCES  FOR  ATTRACTIVE  AND  REPULSIVE  POTENTIALS 
ROEA=3 . 5G/CVR 
ROEB=2oO/CVR 
R0EC2=ROEC*R0EC 


TERS 

N  MAX 

,  AND 

ROEB) 

G. 

A(l 

A(l 

A(l 

A(l 

A(l 

A(2 

A(2 

A(2 

A(2 

A(2 

A(3 

A(3 

A(3 

A(3 

A(3 


FOR  CALCULATION  OF  THE  BEST  CUBIC  FIT  IN  THE  GAP 
IMUM  DISTANCE  CUTOFF  OF  THE  PEPULSIVE  POTENTIAL 
MINIMUM  DISTANCE  CUTOFF  OF  THE  ATTRACTIVE  POTEN- 
SUBPOUTINE  CROSYM  ACTUALLY  PERFORMS  THIS  CURVE 


,1 
,2 
,3 
,4 
,5 
,1 
,2 
,3 
,4 
,5 
,1 
,2 
,3 
1 4 
,5 


)=1.0 

)  =  ROE 
)=ROE 
)=ROE 
)  =  EXP 
)=1.0 
)=ROE 
)=ROE 
)=R0E 
)=EXP 
)=0.0 

)=-2. 


A 

A*ROEA 

A**3 

( EXA+EXB*ROEA) 


B*R0EB 
(CGD1+CGB1*R0EB)-EXP(CGD2+CGB2*R0EB) 


0 

C-ROEA 
0*ROEA*ROEA 
)=EXP(FXA+EXB*ROEA)/CVED 
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c 

SEL 
100 

AND 

IAB 

CRY 

AL 

ANA 


A(4,l)=0.0 

A(4,2)=-1.0 

A(4f3)=-2.0*R0 

A(4»4)=-3.0*R0 

A(4,5)=(EXP(CG 

CALL  CROSYM 

CP0=A(1,5) 

CP1=A(2,5) 

CP2=A(3,5) 

CP3=A<4,5) 

CF0=-CP1*CVEP 

CF1=-2.0*CP2*C 

CF2=-2.0*CP3*C 


EB 

EB*ROEB 

F1+CGB1*RGEB)-EXP(CGF2+CGB2*R0EB)  )/CVED 


ECTION 
i  110, 

DI  AMG 
LES  ES 
STAL. 
X-POSI 
LOGOUS 
GO 


OF  THE  DES 
AND  111  PL 
ND  STRUCTUR 
TAELISHING 

PXI(I)  AND 
TION  GF  THE 


VED 
VED 

IRED  CRYSTAL  STRUCTURE  AND  ORIENTATION. 
ANES  OF  FACE-CENTERED,  BODY-CENTERED, 
ES  ARE  ALLOWED.   ILAY  AND  IDEEP  ARE  VAR- 
THE  NUMBER  OF  MOBILE  LAYERS  IN  THE 
RXK(I)  ARE  VARIABLES  SAVING  THE  ORIGIN- 
I«TH  ATOM.   Y  AND  Z  POSITIONS  ARE 


11 

12 
13 
14 
15 
16 
17 


19 
30 

35 

40 


45 
C 

THIS  S 
DIFFER 
PARAME 
ROUTIN 
INTERS 
TIONS 

50 


CAL 

GO 

CAL 

GO 

CAL 

GO 

CAL 

GO 

CAL 

GO 

CAL 

GO 

CAL 

GO 


TO  (  11 
L  L100 
TO  30 
L  1.110 
TO  30 
L  LI  11 
TO  30 
L  B100 
TO  30 
L  B11C 
TO  30 
L  Bill 
TO  30 
L  D100 
TD  30 

L  Dim 


12, 13, 14, 15, 16, 17, 18, 19), LS 


GO  TO  30 

CALL  Dill 

ILAY=IDEEP 

IF (IDEEP)  3  5,3  5,40 

LD=LL 

ILAY=IY 

DO  45  I=1,LL 

RXK( I )=RX(  I  ) 

RYK( I )  =  RY( I  ) 

RZK( I J=RZ( I ) 

RXI ( 1 )=RX( I ) 

RYI ( I )=RY( I ) 

RZI (I )=RZ(I ) 


55 
60 


ECT ION 
ENT  DAT 
TER  CAL 
E  PLACE 
TITIALS 
IN  THE 
IF(NRUN 
READ  ( 
D1X,D1Y 
IF(DTI. 
DO  55  I 
LCUT( I  ) 
RX( I )=R 
RY( I)=R 
RZ(  I  )  =  P 
RXK( I)= 
RYK( I )= 
RZK( I)= 
NRUN=1 
CALL  PL 
RXI (1)= 


ONE  TO  REPEAT  A  RUN  OF  THE  PROGRAM  WITH 
UT  REPEATING  INITIALIZATION,  POTENTIAL 
NS  AND  CRYSTAL  LATTICE  BUILDING.   SUB- 
CUT(I)  AND  NVAC  TO  CREATE  VACANCIES, 


GO  TO  60 

)  EVR,DTI  ,NTT,NS , ND , I P , I  DEE P , I  TYPE, NVAC, 


ALLOWS 
A  WITHO 
CULATIO 

USES  L 
,  AND  REPLACEMENT  IMPURITIES  AT  DESIRED  LOCA 
LATTICE 
•EQ.t ) 

5,9040 
,D1Z, IQ 

EO.O)  GO  TO  9999 
=  1,LL 
=  0 

XI  (  I  ) 
YI(  I  ) 
Z  I  (  I  ) 
RXI ( I ) 
RYI { I ) 
RZK  I) 

ACE 
RX(1) 
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RYI (1)=RY(1 
RZI(1)=RZ( 1 
RXKU )=RX( 1 
RYK(1)=RY( 1 
RZK(1)=RZ(1 
DO    65    1=1, LL 
VX( I )  =  C.O 
VY( I )=0.0 
VZ( I) =0.0 
PPE( I J=0.0 
PKE( I)=0.0 
65    PTE(I)=0.0 
TPOT=C.O 
NEW=C 
C 

THE  ENERGY  SUBROUTIN 
EACH  ATOM  IN  THE  LAT 
ENERGY  FOR  ALL  ATOMS 
DEFECT. 

CALL  ENERGY 

CALL  LOCAL 

BIND=-TPOT 

TPKE=0.0 

TLPEO=TLPE 

DLPE=TLPE-TLPEi 

TE=TPOT+BIND 


E  CALCULATES  THE  POTENTIAL 
TICE.   SUBROUTINE  LOCAL  SUM 
WITHIN  A  SPECIFIC  RADIUS  0 


ENERGY  OF 
S  UP  THIS 
F  THE  POINT 


C 

THIS  S 
UNITS, 
TIME 


ECTION  PRINTS 
AND  BINDING  E 
ZFRO. 
TIME=0.0 
NT  =  0 

(  6,961 

(  6,962 

1=1, LL, 3 


OUT  X,  Y,  AND  Z  COORDINATES 
NERGIES  OF  EACH  ATOM  IN  THE 


,  IN  LATTICE 
CRYSTAL  AT 


C 

THIS 
FORCE 


WRITE  (  6,9610) 

WRITE  (  6,9620)  IH2,NT 
DO  70 
K=I  +  1 
J  =  1+2 
70   WRJTE  (  6,963 
1RY(K) ,RZ(K) ,PP 

WRITE  (  6,964 
NPAGE=1 
NPAGE=NPAGE+1 

WRITE  (  6,9650)  NPAGE 


G)  I  ,RX(  I  ),RYU  )  ,RZ(  I  )  ,PPE( 
E  (  K  )  ,  J  ,  R  X  (  J  )  ,  R  Y  (  J  )  ,  R  Z  (  J  )  ,  P  P 
C)  TPKE ,TPOT,TE,TLPE,ROEL,D 


IS  THE  MAIN  BODY  OF  THE  PRO 

METHOD,  EXPLAINED  IN  DETAI 
THE  DYNAMICS  FOR  EACH  INDIVIDUAL 
CALCULATES  ALL  MUTUAL  FORCES  AMU 
FORCES,  THIS  SECTION  THEN  CALCUL 
THE  PRIMARY,  AND  ALL  OTHER  ATOMS 
STEP;  AND  THEN  RECALCULATES  FINA 
AND  ALL  OTHER  ATOMS,  BASED  ON  TH 
FORCES.  THIS  SECTION  ALSO  INCLU 
CULATIONS,  BASED  ON  THE  VELOCITI 
CALCULATES  A  NEW  TIMESTEP  DURATI 
STEP,  BASED  ON  EITHER  A  MAXIMUM 
ALLOWED  ENERGY.  (SEE  APP.  C)  V 
END  OF  EACH  TIMESTEP  AS  A  METHOD 
95  TFAC=2.0*PTMAS*DTI*CVD 

TFACB=2 . 0*PBMA S*DT I *C VD 

DT=1.0E-14 
100  DTOD-DT/CVD 

HDT0D=0.5*DT0D 

DTOM=DT/PTMAS 

HDT0H=0.5*DT0M 

DTOMB=DT/PBMAS 

HDT0MB=0.5*0T0 
200  CALL  STEP 

IF(LCUT(  D.GT. 

1=1 

RXK( I )=RX( I ) 

RYK( I )  =  RY< I ) 
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HE  NEXT  TIME- 
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ALVED  AT  THE 


MB 

0)  GO  TO  240 
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RZK( I)=RZ( I ) 

RX(  I  )=RX( I  )+DTOD*(HDTOMB*FX( I  )  +  VX( I  )  ) 

RY( I )=RY( I )+DTOD*(HDTOMB*FY( I)+VY( I ) ) 

RZ( I )=RZ(I )+DTOD*(HDTOMB*FZ< I  )  +  VZ(  I )  ) 
2  40  DO  245  1=2, LO 

IF(LCUT( I ).GT.O )G0  TO  245 

RXK( I J=RX( I ) 

RYK( I  )=RY(  I  ) 

RZK( I )=RZ( I ) 

RX( I)=RX(I  )+DTOD*(HDTOM*FX( I )+VX( I  )  ) 

RY( I)=RY( I )+DTOD*(H0TOM*FY( I )+VY( I ) ) 

RZ( I )=RZ( I )+DTOD*( HDTOM*FZ(  I  )+VZ(  I  )  ) 
245  CONTINUE 

CALL  STEP 

EMAX=0.0 

FMAX=0.0 

TIME=TIME+DT 

NT=NT+1 

IF(LCUT( D.GT.O)  GO  TO  265 

1  =  1 

VSS  =  VX( I  ) 

VX( I)=VSS+HDTOMB*FX(I  ) 

RX( I)=RXK( I )+( VX( I )+VSS)*HDTOD 

VSS=VY(  I) 

VY( I )=VSS+HDTOMB*FY(I ) 

RY( I)=RYK(I )+CVY(I )+VSS)*HDTOD 

VSS=VZ<  I  ) 

VZ( I )=VSS+HDTOMB*FZ(I ) 

RZ( I )=RZK( I )  +  (VZ(I  )+VSS  )*HDTOD 

PKE(  I )=VX<  I )*VX(  I  )+VY(  I)*VY<I)+VZ(I)*VZ(I) 

EMAX=PKE(  I  ) 

FMAX=FX(I )*FX(I)+FY(I)*FY(I)+FZ(I)*FZ(I) 

FORCl=r:MAX 
260  FX(I)=0.0 

FY( I )=C.O 

FZ(  I  1=0.0 
265  Dn  9Q0  I =2 » LD 

IF(LCUT( I).GT.O)GO  TO  2  80 

VSS  =  VX(  I  ) 

VX( I)=VSS+HDTOM*FX(  I  ) 

RX( I )=RXK( I )  +  ( VX(I )+VSS)*HDTOD 

VSS  =  VY(  I  ) 

VY(  I  )-VSS+HDTOM*FY<  I  ) 

RY( I)=RYK( I )+( VY(I )+VSS)*HDTOD 

VSS=VZ(  I  ) 

VZ( I )~VSS+HDTOM*FZ(  I  ) 

RZ( 1)=F  ZK( I )  +  (VZ(I  ) +VSS )*HDTOD 

PKE( I)=VX( I )*VX( I )+VY(  I  )*VY< I ) +VZ (  I  )*VZ ( I  ) 
27  5  F2=FX( I)*FX(I)+FY(I)*FY(I)+FZ(I)*FZ(I) 

FX( I )=0.0 

FY(I)=0.0 

FZ( I >=C.Q 

IF(F2*GT«FMAX)  FMAX=F2 

IF(PKE( I ).GT.EMaX)  EMAX=PKE( I ) 
280  CONTINUE 

IF(EMAXoEQ.O.Q  )  GO  TO  285 

IF(FMAX.EQ.O.O)  GO  TO  285 

GO  TO  287 
285  DT=1.0E-14 

GO  TO  300 
287  IF(EMAX.EQ.PKE( 1  )  )  GO  TO  29U 

IF(FMAX.EQ.FORCl)  GC  TO  290 

EMAX=SQRT( EMAX ) 

FMAX=SQRT( FN!  AX ) 

DT=SQRT(TFAC/FMAX) 

FTERM=FMAX*DT/ ( 2.0*PTMAS) 

GO  TO  295 
290  EMAX=SQRT(PKE( 1 ) ) 

FMAX=SORT(FORC1 ) 

IF(FORCl.EQ.O.C)  GO  TO  285 

DT=SQRT(TFACB/FMAX) 

FTERM=FMAX*OT/ <2.0*PBMAS) 
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295 


300 
310 
3  20 


325 
340 

3  50 

C 

THE  PR 
TION  A 

400 
C 

POTENT 
ARE  CA 
TOTAL 
FORMED 
INITIA 

410 

4  50 


IF(EMAX.LE.FTERM) 
FDTI=DTI*CVD 
DT=FDTI/EMAX 
IF( ISHUT.EQ.-l 
IF(NS-NT)  400, 
DO  325  1=1 ,LL 
VX( I )=0.5*VX( I 
VY(  1  )=0.5*VY(1 
VZ( I )=C.5*VZ( I 
GO  TO  100 
DO  350  1=1,  LL 
RX( I)=RXK( I  ) 
RY(  I)=RYK(  I ) 
RZ(  I  )-P.ZK(  I  ) 


GO  TO  300 


)  GO  TO 
4C0, 320 

I 

) 
) 


400 


620 
700 


C 

THIS  S 
ENERGY 
EVERY 
#NTT. 

720 
1 

7  50 


NTT=NT 

INT  SUBROUTINE  PLACES  A  HEAD 
T  THE  TOP  OF  EACH  TIMESTEP  P 
CALL  PRINT 

IAL  ENERGY  AND  LOCAL  POTENTI 
LCULATED  BASED  ON  THE  NEW  PO 
POTENTIAL  AND  KINETIC  ENERGY 
DX,  DY,  AND  DZ  KEEP  TRACK 
L  POSITION  AT  TIME  ZERO  FOR 
TPOT=0.0 
DO  450  I=1,LL 
PPE( I )=0.0 
PTE( I )=0.0 
CALL  ENERGY 
CALL  LOCAL 
DLPE=TLPE-TLPEO 
PKE( 1)=HBMAS*PKE(1 ) 
TPKE=PKE(1 ) 
PTE(1)=PKE(1)+PPE( 1) 
DO  62  0  1=2, LL 
PKFf  T  )  =  HTMAS*PKF (  I  ) 
TPKE=TPKt+PKE( I ) 
PTE(  I)=PKE(  I  )+PPE(  I ) 
TE=TPOT+BIND 

WRITE  (  6,9660) 
DO  750  1=1, LD 
UX(  I)=RX(  I  )-P.XI  (  I  ) 
DY( I)=RY(  I  )-RYI  ( I ) 
DZ(  I)  =  RZ(  D-RZI  (  I  ) 


ING  OF  PERTINENT 
RINTOUT. 


INFORM- 


AL ENERGY  FOR  EACH  ATOM 
SITIONS.   SUMMATIONS  OF 
FOR  THE  LATTICE  ARE  PER- 
OF  MOTION  RELATIVE  TO  THE 
EACH  ATOM. 


ECTION  PRINTS  THE  RELATIVE  M 

OF  EACH  ATOM,  FOR  EVERY  TIM 

ND'TH  TIMESTEP,  BEGINNING  WI 


OTION,  VELOCITY,  AND 
ESTEP  SO  DESIGNATED:   IE, 
TH  #NS  AND  ENDING  WITH 


WRITE  (  6,967  0)  I , DX (  I  )  , DY (  I  )  , DZ { I ) , VX (  I )  , VY (  I ) , 
VZ( I ) ,PKE( I) ,PPE( I ) ,PTE(I ) 
CONTINUE 

WRITE  (  6,9640)  TPKE,TPOT,T 

WRITE     (     6,9650)     NPAGE 


E,TLPE,ROEL,DLPE 


760 


780 
7  90 

9  50 
C 
THIS    S 

ENERGI 
PROGRA 
955 


NPAGE=NPAGE+1 

IF(NT-NTT)  760,950,950 

DO  780  I=1,LL 

VX(  I)=0.5*VX(I  ) 

VY(  I)=0.5*VY(I  ) 

VZ(  I  )=0.5*VZ(  I  ) 

CONTINUE 

NS=NS+ND 

GO  TO  100 

CONTINUE 

ECTION  PRINTS  OUT  X,  Y,  AND 
ES  OF  EACH  ATOM  IN  THE  CRYST 
M. 

WRITE  (  6,9620)  IH2,NT 
DO  965  1=1, LL, 3 
K=I  +  1 


Z  COORDINATES  AND  BINDING 
AL  AT  THE  END  OF  THE 
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J=I+2 
965   WRITE  (  6,9630)  I , RX ( I ) , RY ( I  )  , RZ (  I  ) , PP E <  I  ) , K ,RX( K ) , 
1RY(K)  ,RZ(K)  ,PPE  (K)  ,  J,RX<  J)  ,RYU)  ,RZ(  J  ),PPE(  J) 
WRITE  (  6,9640)  TPKE , TPOT , Tc , TLPE , RCEL , D LPE 
WRITE  (  6,965C)  NPAGE 
1000  IF(ISHUT)  9999,50,50 
9999  STOP 
END 

SUBROUTINE  CROSYM 
C 

SOLVES  N  SIMULTANEOUS  EQUATIONS  BY  THE  METHOD  OF  CROUT 
THIS  SUBROUTINE  FITS  THE  BEST  CUBIC  BETWEEN  THE  REPULSIVE 
AND  ATTRACTIVE  PARTS  OF  THE  POTENTIAL. 
C 

COMMON/COMA/  A(4,5),MCR0 

M=MCRQ 

N  =  M+1 

11  =  1 
100  13=11 

SUM=ABS(A( 11,11)) 

DO  120  1  =  11, M 

IF(SUM-ABS(A( 1,11)))  110,120,120 
110  13=1 

SUM=ABS(A(  1,11  )  ) 
120  CONTINUE 

IF ( 13-1 1 )  130,150,130 
130  DO  140  J=1,N 

SUM=-A( II, J) 

A(  II, J)=A(  13, J  ) 
140  A( 13, J)=SUM 
150  13=11+1 

DO  160  I=I3,M 
160  A( 1,11 )=A( I, II )/A( 11,11) 
170  ,12=11-1 

13=11+1 

IF(J2)  180,200,180 
180  DO  190  J=I3,N 

DO  190  1=1, J2 
190  A(I1,J)=A(I1,J)-A(I1,I)*A(I,J) 

IF(Il-M)  200,220,200 
200  J2=I1 

11=11+1 

DO  210  I  =  U,M 

DO  210  J=1,J2 
210  A(I,11)=A(I,I1)-A(I,J)*A(J,I1) 

IF(Il-M)  10  0,170,100 
2  20  DO  24  0  1=1, M 

J2=M-I 

I3=J2+1 

A( I3,N)=A( I3,M)/A( 13,13) 

IF(J2)  230,250,230 
230  DO  240  J=1,J2 

2  40  A(J,N)=A(J,N)-A(I3,N)#A(J,I3) 
250  RETURN 

END 

SUBROUTINE  L100 
C 

THIS  IS  A  LATTICE  GENERATOR  FOR  THE  FCC  (ICO)  ORIENTATION, 
THE  CRYSTAL  IS  DEVELOPED  IN  THE  ORDER,  Z  FOLLOWED  BY  Y, 
FOLLOWED  BY  X. 

IT  CONTAINS  A  NONSTANDARD  USE  OF  THE  SURFACE  RELAXATION 
PARAMETER. 
C 

C0MM0N/C0M1/RX( 100C ) ,RY( 1000 ),RZ( 1C00),LCUT( 100C) , 
1L'L,LD,ITYPE,NVAC 

C0MM0N/C0M4/IX,IY, I Z , SC X , SCY , SC Z , IDEEP,D1X,D1Y,D1Z 

DIMENSION  YLAX(2C) 

DATA  YLAX/2C^0.O/ 
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YLAX(l)=-0.097 

YLAX(2)=-0.024 

SCX=1.0 

SCY=1.0 

SCZ=1.0 

M  =  2 

JT  =  0 

Y=-SCY 

DO  60  J=1,IY 

Y=Y+SCY 

KT  =  0 

Z=-SCZ 

DO  59  K=l, IZ 

Z=Z+SCZ 

IT  =  0 

X=-SCX 

DO  58  1=1, IX 

x=x+scx 

ITT=IT+JT+KT 

IF( ITT-(  ITT/2)*2)  57,30,57 

30 

R  X  (  M  )  =  X 

RY(M)=Y+YLAX(J) 

RZ(M)=Z 

M=M+1 

57 

CONTINUE 

IT  =IT+1 

58 

CONTINUE 

KT=KT+1 

59 

CONTINUE 

JT  =  JT  +  1 

IF(  IDFEP-JT)  60,110,60 

60 

CONTINUE 

LL=M-1 

100 

RETURN 

110 

LD=M-1 

GO  TO  60 

CllL' 

SUBROUTINE  L110 
C 

THIS  IS  A  LATTICE  GENERATOR  FOR  THF  FCC  (110)  ORIENTATION. 
THE  CRYSTAL  IS  DEVELOPED  IN  THE  ORDER,  Z  FOLLOWED  BY  Y, 
FOLLOWED  BY  X. 

IT  CONTAINS  A  NONSTANDARD  USE  OF  THE  SURFACE  RELAXATION 
PARAMETER. 
C 

C0MM0N/C0M1/RX( 1000) ,RY( 1000) ,RZ(1000) ,LCUT( 1000) t 
1LL,LD, ITYPE,NVAC 

C0MM0N/C0f-'.4/IX,IY,  I  Z  ,  SC  X  ,  SCY  ,  SCZ  t  IDEEP,D1X,D1Y,D1Z 

DIMENSION  YLAX(20) 

DATA  YLAX/2C*0«0/ 

YLAX( l)=-0.07 

YLAX(2)=-Oo02 

RO=1.0/SQRT(2.0) 

SCX=RO 

SCY=RO 

SCZ=1.0 

M  =  2 

JT  =  0 

Y=-SCY 

DO  60  J=1,IY 

Y=Y+SCY 

KT=0 

Z=-SCZ 

DO  59  K  =  l,  IZ 

Z=Z+SCZ 

IT  =  0 

X=-SCX 

DO  58  1=1, IX 

X= X  +  SC  X 

IF{ IT-( IT/2)*2)  21,11,21 
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11  IF( JT-( JT/2)*2)  57,12,57 

12  IF(KT-(KT/2)*2)  57,30,57 

21  IF< JT-( JT/2)*2)  22,57,22 

22  IF(KT-(KT/2)*2)  30,57,30 
30  RX(M)=X 

RY(M)=Y+YLAX( J) 

RZ(M)=Z 

M  =  M  +  1 

57  CONTINUE 
IT  =  IT+1 

58  CONTINUE 
KT=KT+1 

59  CONTINUE 

JT  =  JT  +  1 
IF(IOEEP-JT)  60,110,60 

60  CONTINUE 
LL=M-1 

100  RETURN 
110  LD=M-1 

GO  TO  60 

END 


SUBROUTINE  Llll 


C 

THIS  IS  A  LATTICE  GENERATOR  FOR  THE  FCC  (111)  ORIENTATION. 

THE  CRYSTAL  IS  DEVELOPED  IN  THE  ORDER,  Z  FOLLOWED  BY  Y, 

FOLLOWED  BY  X. 

IT  CONTAINS  A. NONSTANDARD  USE  OF  THE  SURFACE  RELAXATION 

PARAMETER. 

C 

C0MM0N/CCM1/RX( 1C00),RY( 1G00) , R Z ( 1000 ) , LCUT( 1000) , 
1LL,LD, ITYPE,NVAC 

C0MMCN/CGM4/IX, IY, I Z , SCX , SCY , SCZ , IDEEP,D1X,D1Y,D1Z 

DIMENSION  YLAX(20) 

DATA  YLAX/?0*.T.O/ 

YLAX! l)=-0.04 

YLAX(2)=-0.01 

SCX=1.0/SQRT(2.0) 

SCY=2«0/SQRT(3.0) 

SCZ-SQRTC1.5) 

SSCZ=SCZ/3.C 

M  =  2 

JT=0 

Y=-SCY 

DO  60  J=l, IY 

Y=Y+SCY 

JTS=JT+JT/3 

z=-scz 

KT=0 

DO  59  K=1,IZ 
Z=Z+SCZ 
IT  =  0 
X=-SCX 

00  58  1  =  1,  IX 
X=X+SCX 
1N=IT+JTS+KT 

IF( IN-( IN/ 2)*2 )  57,30,57 
30  RX(M)=X 

RY(M)=Y+YLAX( J ) 

IF( JT-3*(JT/3) )  41,45,41 

41  JTT=JT 

42  JTT=JTT-3 
IF(JTT)  43,45,42 

43  JTT=JTT+3 
ZP=JTT 

RZ(M)=Z+ZP*SSCZ 
GO  TO  50 

45  RZ(M)=Z 
50  M=M+1 
57  CONTINUE 
IT=IT+1 
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58 
59 


60 

100 
110 


CONTINUE 
KT=KT+1 

CONTINUE 

JT=JT+1 

IF( IDEEP-JT) 

CONTINUE 

LL=M-1 

RETURN 

LD=M-1 

GO  TO  60 

END 


60,110,60 


THE  THE  BCC  (100)  ORIENTATION, 
THE  ORDER,  X  FOLLOWED  BY  Z, 


SUBROUTINE  B100 
C 

THIS  IS  A  LATTICE  GENERATOR 
THE  CRYSTAL  IS  DEVELOPED  IN 
FOLLOWED  BY  Y. 

IT  CONTAINS  A  NONSTANDARD  USE  OF  THE  SURFACE  RELAXATION 
PARAMETER. 
C 

COMMON/ COM 1/RX(10C0),RY( 1000 ),RZ( 1000),LCUT( 1000) , 
1LL,LD, ITYPE,NVAC 

C0MM0N/CCM4/IX,IY. I Z , SCX , SCY , SCZ , I  DEEP , Dl X  ,  Dl Y, Dl Z 

DIMENSION  YLAX (20) 

DATA  YLAX/2C'*0.0/ 

YLAX( 1)=-0.2G 

YLAX(2)=-0.03 

SCX=1.0 

SCY=loO 

SCZ=1.0 

M  =  2 

JT  =  0 

Y=-SCY 

DO  60  J=1,IY 

Y=Y+SCY 

i/T-ri 

z=-scz 

DO  59  K=1,IZ 

Z=Z+SCZ 

IT=0 

X=-SCX 

DO  58  1  =  1,  IX 

X  =  X+  SC  X 

IF( lT-(  IT/2)*2)  21,11,21 

11  IF( JT-( JT/2)*2)  57,12,57 

12  IF(KT-(KT/2)*2)  57,30,57 

21  IF(JT-( JT/2)*2)  22,57,22 

22  IF(KT-(KT/2)*2)  30,57,30 
30  RX(M)=X 

RY(M)=Y+YLAX( J) 

RZ(M)=Z 

M  =  M  +  1 

57  IT=IT+1 

58  CONTINUE 
KT=KT+1 

59  CONTINUE 

JT  =  JT  +  1 
IF(IDEEP-JT)  60,110,60 

60  CONTINUE 
LL=M-1 

100  RETURN 
110  LD=M-1 

GO  TO  60 

END 
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SUBROUTINE  BllO 

IS  A  LATTICE  GENERATOR 

RYSTAL  IS  DEVELOPED  IN 

WED  BY  Y. 

NTAINS  A  NONSTANDARD  USE 

ETER. 


THE  THE  BCC  (11C)  ORIENTATION, 
THE  ORDER,  X  FOLLOWED  BY  Z, 

OF  THE  SURFACE  RELAXATION 


COMMON/CCM1/RXC1000J t RY( 1000 ) , RZ ( 100C  ),LCUT(  10C0) , 
1LL,LD, ITYPE,NVAC 
COMMON/ C0M4/ IX, IY, I Z, SCX , SCY , SCZ , IDEEP,D1X,D1Y,D1Z 
DIMENSION  YLAX(20) 
DATA  YLAX/2G*C.O/ 
YLAX(1)=-0.10 
YLAX(2)=-0.C1 
SCX  =  SQRT(2.C) 
SCY  =  SQRT(2cO) 
SCZ=1.C 
M=2 
JT=0 
Y=-SCY 

DO  60  J  =  l,  IY 
Y=Y+SCY 
KT=0 
Z=-SCZ 

DO  59  K=l, IZ 
Z=Z+SCZ 
IT=0 
X=-SCX 

DO  58  1=1, IX 
X=X+SCX 
ITT=IT+JT+KT 

IF< ITT-(  ITT/2)*2)     57,30,57 
^P)     RX(M)=X 

DVIM)rV+YI    AX( J ) 

RZ(M)=Z 

M    =    M+l 

57  CONTINUE 
IT  =IT+3 

58  CONTINUE 
KT=KT+1 

59  CONTINUE 

JT  =  JT  +  1 
IF(IDEEP-JT)  60,110,60 

60  CONTINUE 
LL=M-1 
RETURN 

110  LD  =  M-l 
GO  TO  60 
END 


SUBROUTINE  Bill 
C 

THIS  IS  A  LATTICE  GENERATOR  FOR  THE  BCC  (111)  ORIENTATION. 
THE  CRYSTAL  IS  DEVELOPED  IN  THE  ORDER,  X  FOLLOWED  BY  Z, 
FOLLOWED  BY  Y. 

IT  CONTAINS  A  NONSTANDARD  USE  OF  THE  SURFACE  RELAXATION 
PARAMETER. 
C 

C0MM0N/CCM1/RX( 1G0C ) ,RY( 1000 ) ,RZ( 1000 ) ,LCUT( 1000) , 
1LL,LD, ITYPE ,NVAC 

COMMON/ C0M4/ IX, IY, I Z , SCX , SCY , SCZ , IDEEP,D1X,D1Y,D1Z 

DIMENSION  YLAXC20) 

DATA  YLAX/2C*0.0/ 

YLAXQ  )=-0.10 

YLAX<2)=-0.01 

SCX  =  SQRT(2.0) 
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SCY  =  SQRTd. 0/3.0) 

SCZ  =  SQRK6.0) 

M  =  2 

JT=0 

Y=-SCY 

DO  60  J=1,IY 

Y=Y+SCY 

KT  =  0 

JJ  =  J-( JT/3)*3 

GO  TO  (5,10  ,15 ) , JJ 
5  Z  =-SCZ 

GO  TO  20 
10  Z  =-SCZ+SCZ/3.0 

GO  TO  20 
15  Z  =-SCZ+(2.0*SCZ/3.0) 
20  DO  59  K  =  l,  IZ 

Z=Z+SCZ 

IT  =  0 

X=-SCX 

DO  58  1=1, IX 

X=X+SCX 

ITT=IT+JJ+KT-1 

IF(  ITT-(  ITT/2)*2)     57,30,57 
30    RX(M)=X 

RY(M)=Y+YLAX< J) 

RZ(M)=Z 

M    =    M+l 

57  IT=IT+1 

58  CONTINUE 
KT=KT+1 

59  CONTINUE 

JT  =  JT  +  1 
IF(IDEEP-JT)  60,110,60 

60  CONTINUE 
LL=M-1 
RETURN 

110  LD=M-1 

GO  TO  60 
END 

SUBROUTINE  D100 

C0MM0N/C0M1/RX( 10001 ,RY(1000) ,RZ(1000) , LCUT(IOCO) , 
1LL,LD, ITYPE,NVAC 
CCMM0N/CCM4/IX, I Y, I Z , SCX , SCY , SCZ , IDEEP,D1X,D1Y,D1Z 
RETURN 
END 

SUBROUTINE  DUO 

C0MM0N/C0M1/RX( 1000),RY( 1000 ), RZ ( 1000 ), LCUT( 1000) , 
1LL,LD, I  TYPE, NV AC 
C0MM0N/C0M4/IX,IY,IZ,SCX,SCY,SCZ, I  DEEP , Dl X , Dl Y, Dl Z 
RETURN 
END 

SUBROUTINE  Dill 

COMMON/ C0M1/RX( 1000 ), P Y ( 1C00 ), RZ ( 1000 ), LCUT ( 1000) , 
1LL,LD,ITYPE,NVAC 
COMMON/ C0M4/ IX, IY, I Z , SCX , SCY , SCZ , I  DEE P , Dl X , Dl Y, Dl Z 
RETURN 
END 

SUBROUTINE  PLACE 
C 

THIS  SUBROUTINE  LOCATES  A  VACANCY,  INTERSTITIAL,  OR  REPLACE- 
MENT IMPURITY  IN  THE  LATTICE. 
C 

COMMON /C0M1/RX( lOOu) ,RY( 1000) ,RZ ( 1G0C ) , LCUT( 1000) , 
1LL,LD, ITYPE,NVAC 
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COMMON/C0M4/IX ,  IY, I Z,SCX, SCY,SCZ, IDEEP,D1X,D1Y,D1Z 

GO  TO   ( 10,20,30),  ITYPE 
10  LCUT(NVAC)  =  1 

LCUT(l)  =  1 

RX(1)=0.0 

RY(1)=0.0 

RZ( 1)=C0 

GO  TO  40 
20  RX(1)  =  RX(NVAC)  -  D1X 

RY(1)  =  RY(NVAG)  -  D1Y 

RZ( 1 »  =  RZ(NVAC)  +  D1Z 

GO  TO  40 
30  LCUT(NVAC)  =  1 

RX(1)  =  RX(NVAC) 

RY(1)  =  RY(NVAC) 

RZ(1)  =  RZ(NVAC) 
40  RETURN 

END 

SUBROUTINE    STEP 
C 

THIS    SUBROUTINE    DOES     THE     DYNAMICS    FOR    ONE    TIMESTEP. 
THE    FIRST    HALF    DOES    THE    DYNAMICS    FOR     ATOM    #1;     THE    SECOND 
HALF    FOR    ALL    OTHERS. 
C 

COMMON/CCMl/RXdGOO)  ,RY(  1000)  ,RZ(  1000)  ,LCUT(  1000)  , 
ILL, LD, ITYPE, NV AC 

COMMO\!/CCM5/ROE,ROE2,ROEM,EXA,EXB,PEXA,PEXB,FXA,PFXA, 
1IQ,TSAVE,BSAVE 

COMMOi\!/COM6/FX(1000)  ,FY(100U)  ,FZ(  1000)  ,  PAC  ,  PFPTC  ,  FM 

COMMON/COM8/ROEA,ROEB,ROEC,RUEC2,CPO,CP1,CP2,CP3, 
1CF0.CF1 ,CF2,CGD1,CGD2,CGB1,CGB2,CGF1,CGF2 

IF( IQ.EQ.l)    GO    TO    200 

1  =  1 

IF  (  LCUT  i  I  »  )  ?00  .1  05,200 
105  IP-1-*-! 

DO  195  J=IP,LL 

IF(LCUTU))     195,110,195 
110    DRX  =  RX( J )-RX(I  ) 

IF(DRX)     113,117,117 
113     IF<DRX+ROE)     195,195,120 
117     IF(DRX-ROE)     12J, 195,195 
120    DRY=RY( J)-RY(I ) 

IF(DP.Y)     123,127,127 
123     IF(DRY+ROE)     195,195,130 
127     IF(DRY-ROE)     130,195,195 
130    DRZ=RZ( J)-RZ(I  ) 

IF(DRZ)     133,137,137 
133    IF(DRZ+ROE)     195,195,140 
137     IF(DRZ-ROE)     140,195,195 
140     DIST=DRX*DRX+DRY*DRY+DRZ*DRZ 

IF(DIST-R0E2)     150,195,195 
150    DIST=SQRT(DIST) 
160     IF(DIST-POEM)     162,162,165 
162    FGRCE=EXP(PFXA+PEXB*DIST) 

GO  TO  180 
165  DFF=ROE-DIST 

IF(DFF-l.OE-lO)     195,195,167 
167    FORCE=(EXP(PAC+PEXB*DIST)-PFPTC)/DFF 
180    IF(FM-FORCE)     190,190,195 
190    FCD=FORCE/DIST 

FA=FOD*DRX 

FX( J)=FX( J)+FA 

FX(  I  )  =  FX(I  )-FA 

FA=FOD*DRY 

FY( J)=FY( J)+FA 

FY(  I )  =  FY( I  )-FA 

FA=FOD*DRZ 

FZ( J  )=FZ( J)+FA 

FZ(  I  )  =  FZ{ I )-FA 
195    CONTINUE 
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200  DO  300  I=IQ,LD 

IF(LCUT( I ) )  300,205,300 
205  IP=I+1 

DO  295  J=IPtLL 

IF(LCUT( J) )  295,210,295 
210  DP.X=RX(  J)-RX(  I  ) 

IF(DRX)  213,217,217 
213  IF(DRX+RCEC)  295,295,220 
217  IF(DRX-ROEC)  220,295,295 
220  DRY=RY( J)-RY( I ) 

IF(DRY)  223,227,227 
223  IF(DRY+ROEC)  295,295,230 
227  IF(DRY-ROEC)  23Ct295,295 
230  DRZ=RZ< J)-RZ(I ) 

IF(DRZ)  233,237,237 
233  IF(DRZ+ROEC)  295,295,240 
237  IF(DRZ-ROEC)  240,295,295 
2  40  D1ST=DRX*DRX+DRY*DRY+DRZ*DRZ 

IF(DIST-R0EC2)  250,295,295 
250  DIST=SORT(DIST) 

IF(DIST-ROEA)  260,255,255 
255  IF(DIST-ROEB)  265,27u,270 
260  FORCE=EXP(FXA+EXB*DIST) 

GO  TO  2  30 
2  65  F0RCE=DIST*(DIST*CF2+CF1)+CF0 

GO  TO  280 
2  70  F0RCE=FXP(CGF1+CGB1*DIST)-EXP(CGF2+CGB2*DIST) 
280  IF(ABS( FORCE). Lb. FM)  GO  TO  295 

FOD  =  FORCE/DIST 

FA=FOO*DRX 

FX( J)=FX( J)+FA 

FX( I )  =  FX( I  )-FA 

FA=FOO*DRY 

FY( J)=FY( J)+FA 

FY( I )=FY( 1 )-FA 

FA=F0D*nR7 

FZ(J  )=FZ( J)+FA 

FZ( I )=FZ( I )-FA 
295  CONTINUE 
300  CONTINUE 

RETURN 

END 

SUBROUTINE  ENERGY 
C 

THIS  SUBROUTINE  CALCULATES  THE  MUTUAL  POTENTIAL  ENERGIES. 
THE  FIRST  HALF  DOES  THE  DYNAMICS  FOR  ATOM  #1;  THE  SECOND 
HALF  FOR  ALL  OTHERS. 
C 

COMMON/ C0M1/RX( 10C0) , PY(IOOO) ,RZ( 1000) , LCUT( 10CC) i 
1LL,LD,  I  TYPE, NV AC 

COMMON /COM5/ROE,ROE2,R0EM,EXA,EXB,PEXA,PEXB,FXA,PFXA, 
1IQ,TSAVE, BSAVE 
COMMON/ C0M7/PPTC,TPOT,PPE( 1000 ) , TLPE , ROEL , R0EL2 , NEW 
C0MM0N/C0M8/R0EA,R0EB,R0EC,R0EC2,CP0,CP1 ,CP2,CP3, 
1CF0,CF1,CF2,CGD1,CGD2,CGB1,CGB2,CGF1,CGF2 
IF(IQ.EQ.l)  GO  TO  2G0 
1  =  1 

IF(LCUT( I) )  600,505,600 
505  IP=I+1 

DO  595  J=IP,LL 
IF(LCUTU))  595,510,595 
510  DRX=RX( J)-RX( I ) 

IF(DRX)  513,517,517 
513  IF(DRX+ROE)  595,595,520 
517  IF(DRX-ROE)  520,595,595 
520  DRY=RY( J )-RY(I ) 

IF(DRY)  523,527,527 
523  IF(DRY+ROE)  595,595,530 
527  IF(DRY-ROE)  530,595,595 
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530  DRZ=RZ( J)-RZ(I ) 

IF(DRZ)  533,537,537 
533  IF(DRZ+R0E)  595,595,540 
537  IF(DRZ-ROE)  540,595,595 
540  DIST=DRX*DRX+DRY*DRY+DRZ*DRZ 

IF(DIST-R0E2)  550,595,595 
550  DIST=SQRT(DIST) 
560  POT=EXP(PEXA+PEXB*DIST)-PPTC 
580  TPOT=TPOT+POT 

PPE(I)=PPE(I )+BSAVE*POT 

PPE( J)=PPE( J)+TSAVE*POT 
595  CONTINUE 
600  CONTINUE 
C 

200  DO  300  I=IQ,LL 

IF(LCUT(  I)  )  300,205,300 
205  IP=I+1 

DO  295  J=IP,LL 

IF(LCUT(J))  295,210,295 
210  DRX=RX( J)-RX( I ) 

IF(DRX)  213,217,217 
213  IF(DRX+ROEC)  295,295,220 
217  IF(DRX-ROEC)  220,295,295 
2  20  DRY=RY( J)-RY( I ) 

IF(DRY)  223,227,227 
223  IF(DRY+ROEC)  295,295,230 
227  IF(DRY-ROEC)  230,295,295 
230  DRZ  =  RZ( J)-RZ( I  ) 

IF(DRZ)  233,237,237 
233  IF(DRZ+ROEC)  295,295,240 
237  IF(DRZ-RCEC)  240,295,295 
240  DIST  =  ORX*DRX+DRY-DRY-f  DRZ*DRZ 

IF(DIST-R0EC2)  250,295,295 
250  DIST=SQRT(DIST) 

IF(DIST-ROEA)  260,255,255 
255  TF(niST-RnFB)  265,270, 27U 
260  POT=EXP( F*A+EYB*nT ST) 

GO  TO  280 
2  65  P0T=DIST*(DIST*(DIST*CP3+CP2)+CP1)+CP0 

GO    TO    2  80 
270    P0T=EXP(CGD1+CGB1*DIST)-EXP(CGD2+CGB2*DIST) 
280    TP0T=TP0T+P0T 

SAVE=0.5*P0T 

PPE(  I)=PPE(  D+SAVE 

PPE( J) =PPE( JJ+SAVE 
295  CONTINUE 
300  CONTINUE 

RETURN 

END 

SUBROUTINE  LOCAL 
C 

THIS  SUBROUTINE  CALCULATES  THE  TOTAL  POTENTIAL  ENERGY  IN  A 
SMALL  VOLUME  AROUND  A  VACANCY  OR  INTERSTITIAL. 
C 

COMMON/ CCM1/PXQ00C  )  ,RY(  1000  )  ,RZ(  1000)  ,LCUT(  1000)  , 
1LL.LD, ITYPE,NVAC 
CCMM0N/CCM7/PPTC,TP0T, PPE( 1000) , TLPE , ROEL , R0EL2, NEW 
D I  MENS  I  ON  LCCAT(50C  ) 
K=l 

TLPE=0.0 

IF(NEW.EQ.1)G0  TO  305 
8  GO  TO  ( 10,20,20) , ITYPE 
10  I=NVAC 

GO  TO  200 
20  1  =  1 
200  DO  300  J=1,LD 

1FU.EQ.J)  GO  TO  250 
210  DRX  =  RX( J  )-RX( I ) 

IF(DRX)  213,217,217 
213  IF(DRX+RGEL)  295,295,220 
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217 
220 

223 
2  27 
2  30 

233 
237 
240 


2  50 

295 
300 

305 

310 


IF(DRX- 

DRY=RY( 

IF(DRY) 

IF(DRY+ 

IF(DRY- 

DRZ=RZ( 

IF(DRZ) 

IF(DRZ+ 

IF(DRZ- 

DIST=DR 

R0EL2=R 

IF(DIST 

L0CAT(  K 

K  =  K  +  1 

CONTINU 

CONTIiMU 

KF=K-1 

DO  31G 

TLPE=TL 

CONTINU 

NEW  =  1 

RETURN 

END 


ROEL 

J)-R 

223 

ROEL 
ROEL 
J  )-R 
233 
ROEL 
ROEL 
X-DR 
OEL* 
-ROE 
)=J 

E 
E 


K=1,KF 

PE+PPE( LOCAT(K) 

E 


)  220,295,295 

Y(  I  ) 

,227,227 

)  295,295,230 

)  230,295,295 

Z(  I  ) 

,237,237 

)  295,295,240 

)  240,295,295 

X+DRY*DRY+DRZ*DRZ 

ROEL 

L2)  250,295,295 


C 

THIS 
TION 
C 


SUBROUTINE  PRINT 


SUBROUTINE 
AT  THE  TOP 


PRINTS  TH 
OF  EACH  T 


E  HEADING  OF  ALL  PERTINENT  INFORMA- 
IMESTEP  PRINTOUT. 


9710 
9720 

9730 
9740 
9741 

9742 

9750 

97  60 
9765 
9770 

9780 


COMMON/ 
1LL,LD,I 

COMMON/ 
1TARGET( 

COMMON/ 
] NT, T IMF 

CO NIMQfvl  / 

COMMON/ 
1  IQ,TSAV 

COMMON/ 
1CFCCF1 

FORMAT ( 

FORMAT ( 
1  UNIT  = 

FORMAT ( 
IE  TEMP 
1V/) 

FORMAT ( 
1  F5.2,2 
1  ), ,  4X 

FORMAT ( 
1  F5.2,2 
1  )  ,  ,  4X 
13.2H)  FR 

FORMAT  ( 
1  F5.2,2 
1  ),  ,  4X 

FORMAT ( 
1F5.2,5H 
110X,4HI 

FORMAT ( 
1F9.5,2X 

FORMAT ( 
1A    =  ,F9. 

FORMAT ( 
1TIAL  PA 
1F10.3, ' 
1E10.3, ■ 

FORMAT ( 
1SE  POTE 
1F8.4, « - 


C0M1 
TYPE 
COM2 
4),T 
COM3 
•  DT. 
rr.M/x 

COM  5 
E,BS 
CCM8 
,CF2 
40X, 
OH  T 
,F7. 
4X,6 
=  F5. 

2H  ( 
1HKE 
,  16 
2H  ( 
1HKE 
,  15 
CM  S 
2H  ( 
1HKE 
i  20 
30  H 
i  Z 
Q  =1 
12H 
,5HP 
12X, 
5/) 
■  WH 
RAME 
i  CP 
t  CF 

«  cu 

NTIA 
CGD 


/RX( 

,NVA 

/IH1 

MAS  , 

/RXI 

DTI  , 

/IX, 

/ROE 

AVE 

/ROE 

,CGD 

10A4 

ARGE 

4,4H 

HMAS 

2,7H 


1000),RY(1000),RZ( 1000 ),L CUT ( 1000) , 

C 

(20) 

BULL 

(100 

TLAY 

ty.  I 

,ROE 


,IH2(8) ,IHS( 10) , IHB( 

ET(4) ,BMAS, PLANEtTEM 

(0) ,RYI(100C ) ,RZI(100 

v.  <;ry  Arv,qr7.  jdffp- 
:2iR0EMtEXA^EXB,PEXAi 


6) , IHT(6) , 

P, THERM 

0 ) ,CVR,EVR, 


ni y . Div. ni 7 
PEXBjFXA^PFXA, 


1,CP2,CP3, 
,F2 


A,ROEB, ROEC,ROEC2,CP0,CP 

1,CGD2,CGB1,CGB2,CGF1,CG 

,/,28X,2CA4, /) 

T    -,4A4,1CHPRIMARY    -     ,4A4 , 1 X , 1 4HL ATT  I CE 

ANG) 
S    =,F7.2,13X,6HMASS    =,F7 

DEG    K, ,18H       THERMAL    CUT 


«2,9X, 14HLATTIC 
OFF    =,F5e2,3H    F 


,A4,8H) 
V,       CRYS 
HVACANCY 
,A4,8H) 
V,       CRYS 
HINTERST 
ITE     ,14/ 
,A4, 8H) 
V,        CRYS 
HRE PLACE 
PRIMARY 
=  ,F5.2,5 
12/  ) 

POTENTIA 
FXA=,F9. 
6A4,3X,  5 

EN«  ,F8.4 
TERS  ARE 
2    =•  ,F10 

1  =  •  ,  E 1 0 
T-OFF  AT 
L    PARAME 

2  =•  ,F8. 


PLANE, , 18H 
TAL  SIZE  ( 
IN  SITE  , 
PLANE, ,18H 
TAL  SIZE  ( 


PRIMARY 
,  12, 3H  X 
14/) 

PRIMARY 


ENERGY  =, 

,  12, 3H  X  ,  12 ,3H 


ENERGY  =, 
,I2,3H  X  ,12  ,3H 
(-,F5,2,2H,-,F5c2,2H,+,F5.2, 


ITIAL 

) 

PLANE, ,18H   PRIMARY 

TAL  SIZE  (  ,I2,3H  X 

MENT  IN  SITE  ,  14/) 

START  POINT  (LU)   X 

X,I3,«  LAYERS  ARE  FR 

L   ,6A4,3X, 5HPEXA=,F 

5) 

HEXA  =,F9.5,2X,5HEXB 


<  R  <• 


// 


.3, 


CGF1  =■ ,F8.4, « 


,F8,4,«  THE 
CPO  =', 

CP3  ='  ,F10.3,/ 
.3,  •  ,  CF2  ='  ,E10.3,/ 
1 , F5.2, • ,  WHEN  R  >  • 
TERS  ARE* ,  8A4,//,10 
4.  «  7  CGB1  =«  ,c8a4, '  , 
CGF2  =• ,F8.4, //) 


ENERGY  =, 

, 12, 3H  X  , 12, 3H 

=,F5.2,5H,  Y  =, 
EE  TO  MOVE ' , 

9. 5,2X,5HPEXB=, 

=,F9.5,2X,5HFX 

MATCHING  POTEN 
F10.3, • ,  CP1  =' 
,  «  CFG  =' 
/) 

,F6,3,'  LU,  MOR 
X,"    CGD1  =', 

CGB2  ='•   ,  F8.4, 
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9790  FORMAT (10H  TIMESTEP  , 1 4, 22 X , 6HDT I  =  ,  F5.3,  5H  LU,  , 

122H   ELAPSED  TIME  (SEC)  =  ,  E10.4,21H,  LAST  TIMESTEP  WA 
IS  =  ,  E10.4/) 

WRITE  (  6,971C)  IHStlHl 

WRITE  (  6,9720)  TARGET, BULL ET, CVR 

WRITE  (  6,9730)  TMAS, BMAS , TEMP , THERM 
GO  TO  (401,402,403),  I  TYPE 

401  WRITE  (  6,9740)  PL ANE , EVR, I  X, I Y, I Z , NVAC 
GO  TO  405 

402  WRITE  (  6,9741)  PL ANE ,E VR , I  X , I Y , I Z , D1X , 01 Y , Dl Z, NVAC 
GO  TO  405 

403  WRITE  (  6,^742)  PL ANE , EVR , I  X , I Y , I Z , NVAC 

405   WRITE  (  6,9750)  R X  I ( 1 )  , RY I ( 1 )  ,RZ I ( 1 )  ,  I  LAY  ,  I Q 

WRITE  (  6,9760)  IHB, PEXA, PEXB , PFXA 

WRITE  (  6,9765)  I HT , EXA , EXB , FXA 

WRITE  (  6,9770)  ROEA, ROEB , CPO ,CP1 , CP2, CP3 ,CF0 ,CF 1 , CF2 

WRITE  (  6,9780)  ROEC , ROEB , I H2 , CG01 , CGD2 , C GB1 , CGB2  , 

WRITE  (  6,9790)  NT  ,DT I , T IME , DT 
RETURN 
END 
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